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ABSTRACT 

We present global, three-dimensional numerical simulations of HD 189733b and HD 209458b that couple 
the atmospheric dynamics to a realistic representation of non-gray cloud-free radiative transfer. The model, 
which we call the Substellar and Planetary Atmospheric Radiation and Circulation (SPARC) model, adopts the 
MITgcm for the dynamics and uses the radiative model of McKay, Marley, Fortney, and collaborators for the 
radiation. Like earlier work with simplified forcing, our simulations develop a broad eastward equatorial jet, 
mean westward flow at higher latitudes, and substantial flow over the poles at low pressure. For HD 189733b, 
our simulations without TiO and VO opacity can explain the broad features of the observed 8 and 24-/^m light 
curves, including the modest day-night flux variation and the fact that the planet/star flux ratio peaks before the 
secondary eclipse. Our simulations also provide reasonable matches to the Spitzer secondary-eclipse depths 
at 4.5, 5.8, 8, 16, and 24 /im and the groundbased upper limit at 2.2 fim. However, we substantially underpre- 
dict the 3.6 fim secondary-eclipse depth, suggesting that our simulations are too cold in the 0.1-1 bar region. 
Predicted temporal variability in secondary-eclipse depths is ^ 1 % at Spitzer bandpasses, consistent with re- 
cent observational upper limits at 8 ^m. We also show that nonsynchronous rotation can significantly alter the 
jet structure. For HD 209458b, we include TiO and VO opacity; these simulations develop a hot (> 2000 K) 
dayside stratosphere whose horizontal dimensions are small at depth but widen with altitude. Despite this 
stratosphere, we do not reproduce current Spitzer photometry of this planet. Light curves in Spitzer bandpasses 
show modest phase variation and satisfy the observational upper limit on day-night phase variation at 8 /im. 
Subject headings: planets and satellites: general, planets and satellites: individual: HD 209458b, methods: 
numerical, atmospheric effects 



1. INTRODUCTION 

Blasted by starlight 10^-10^ times stronger than that re- 
ceived by Jupiter and experiencin g modest rot ation rates 
due to their presumed tidal locking dGuillot et al.lil99 6), hot 
Jupiters occupy a fascinating meteorological regime that does 
not exist in our Solar System (for an extensive review see 
IShowman et al.l l2008bh . Despite the wide range of transit- 
ing exoplanets that have been discovered, HD 189733b and 
HD 209458b remain the best characterized hot Jupiters and 
represent important test cases for our understanding of these 
objects generally. A variety of observations now exist that 
constrain the 3D temperature structure, composition, hazes, 
and albedo for these two planets. There is now hope that, by 
comparing these observations with detailed models, insights 
into this novel atmospheric regime can be achieved. 

The strongest observational evidence for atmospheric mo- 
tions comes from infrared light curves obtained with the 
Spitzer Space Telescope. Continuous light curves of HD 
189733b over half an orbit have now been obtained at both 8 
and 24 nm, constrairiing th is planet's day-night heat transport 
dKnutson et al.ll2007l l2009i) . The inferred dayside and night- 
side brightness temperatures are ~ 1250K and ^ lOOOK, re- 
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spectively.^ Because the nightside temperatures would be ex- 
tremely cold 200 K) in the absence of winds, these obser- 
vations imply the existence of a vigorous atmospheric circula- 
tion that efficiently transports heat from dayside to nightside. 
Further evidence for winds is the fact that the peak flux in both 
light curves occurs before secondary eclipse, implying that 
the hottes t region hes 20 -30° of long itude east of the substel- 
lar point (iKnutson et al.ll2007l l2009b . For HD 209458b, cur- 
rent data contain insufficient temporal sampling to determine 
whether similar offsets exist but nevertheless demonstrate that 
the 8-/i m day-night brightn ess temperature difference is also 
modest jCowan et al.ll2007l) . 

For both planets we now also have dayside photometry 
at all Spitzer broadband channels (centered at 3.6, 4.5, 5.8, 
8.0, 16, and 24 /im), placing constraints on dayside compo- 
sition and vertical temperature profile. These data were ob- 
tained by differencing the photometry of star+planet taken 
just before/after secondary eclipse from photometry taken 
during secondary eclipse, when only the star is visible. When 
compared with one-dimensional (ID) atmosphere models, 
these data suggest that, near the photosphere pressures, the 
daysi de temperature of HD 1 89733b decreases with alti- 
tude dCharbonneau et al.ll2008t lBarmanll2008t iKnutson et all 
l2009l) . whereas HD 2094 58b instead contains a thermal inver- 
sion ( a hot stratosphere) (IKnutson et al.ll2008t iBurrows et al.l 
l2007h . 

Following pione erin g work by iHubeny et al.l (l2003l) . 
iFortney etal](l2008h and lBurrows et al.l(l2008h suggested that 



* At 8/im, the dayside and nightside brightness temperatures are 1258 ± 
1 1 K and 101 1 ± 5 1 K. At 24 /tm, the dayside and nightside brightness tem- 
peratures are 1220 ± 47 K and 984 ± 48 K. 
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hot Jupiters subdivide into two classes depending on whether 
or not their atmospheres contain highly absorbing substances 
such as gaseous TiO and VO. For a Sun-like primary, solar- 
composition planetary atmospheres inward of 0.04-0.05 AU 
are hot enough to contain TiO and VO; because of the extreme 
visible-wavelength opacity of these compounds, such planets 
absorb starlight at low pressure (^mbar) and naturally exhibit 
dayside stratospheres. For planets outward of ~ 0.05 AU, 
temperatures drop sufficiently for TiO and VO to condense in 
the deep atmospheres; these planets absorb starlight deeper in 
the atmosphere, lack stratospheres, and show spectral bands in 
absorption. Based on simple cornpariso ns of radiative and ad- 
vective timescales. lFortney et al.l(l2008h further suggested that 
the former category (dubbed "pM" class planets) would ex- 
hibit large day-night temperature differences whereas the lat- 
ter category ("pL" class) would exhibit more modest temper- 
ature contrasts. Their calculations suggest that HD 189733b 
is a pL-class planet while HD 209458b is a pM-class planet. 
This dichotomy makes these two planets a particularly inter- 
esting pair for comparison. 

Given these developments, there is a pressing need for 3D 
calculations of the atmospheric circulation on hot Jupiters. 
Several groups have investigated a range of 2D and 3D 
models S howman & Guillot 2002t ICho et all 12 003. 2008; 



Cooper & Showman 2005, 2006; Showman etal. 2008a b; 
Langton & Laughlin 2007, 2008; Dobbs-Dixon & Lin 2008). 
However, to date, all published models have adopted severe 
approximations to the radiative transfer or excluded radiative 
heating/cooling entirely. While such simplified approaches 
are invaluable for investigating the underlying dynamics, a 
detailed attempt to explain wavelength-dependent photome- 
try and light curves must include a realistic coupling of the 
radiative transfer to the dynamics. 

Here, we present new numerical simulations that couple a 
realistic representation of non-gray cloud-free radiative trans- 
fer to the dynamics. Surprisingly, coupling of 3D dynamics 
to nongray radiative transfer has never previously been done 
for any giant planet, not even Jupiter, Saturn, Uranus, and 
Neptune. This is the first 3D dynamical model for any giant 
planet — inside our Solar System or out — to do so. We dub 
our model the Substellar and Planetary Atmospheric Radia- 
tion and Circulation model, or SPARC model, and to honor its 
dynamical heritage usually refer to it as SPARC/MITgcm. §2 
presents the model. §3 describes the basic circulation regime 
and resulting light curves and spectra for HD 189733b. §4 de- 
scribes the effect of a nonsynchronous rotation rate, §5 con- 
siders HD 209458b, and §6 concludes. 

2. MODEL 

In the absence of dynamics, radiative transfer drives tem- 
peratures toward radiative equilibrium, which for a tidally 
locked irradiated planet means hot on the dayside and cold 
on the nightside. These thermal contrasts produce horizon- 
tal pressure gradients that generate winds, which drive the 
atmosphere away from radiative equilibrium. Because of 
this dynamical response, the radiative heating/cooling rate 
is nonzero, and it is this heating/cooling rate that drives the 
dynamics. In our previous work, we calculated the heat- 
ing rate using a simplified Newtonian relaxation scheme 
jShowman & Guillot"2002 HCooper & Showmanll200l 120061: 
[showman et al. 2008a), but here we instead self-consistently 
calculate the heating rate from the radiative transfer. Fig- 
ure[T]illustrates the coupling between dynamics and radiation 
in SPARC/MITgcm. 




Dynamics 

Solves for new v, T 
given heating rate 
and old V, I 



Radiative Transfer 

Solves for upwaid and 
dowmvard fluxes given 
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Fig. 1 . — Coupled interrelationship between dynamics and radiation in our 
general circulation model (GCM), the SPARC/MITgcm. Spectea are calcu- 
lated offline from GCM output. 

2.1. Dynamics 

We solve the global, 3D p rimitive equations in spherical 
geometry using the MITgcm dAdcroft et al.ll2004l) . which is 
a state-of-the-art atmosphere and ocean circulation model 
maintained at MIT. The primitive equations are the stan- 
dard equations for atmospheric flows in statically stable at- 
mospheres when the horizontal dimensions greatly exceed 
the vertical dimensions. For HD 189733b and HD 209458b, 
with horizontal dimensions of 10^-10^ m and scale heights of 
^ 200km a nd 500 km, respectively, we expect aspect ratios of 
~ 20-500. Showman et al . (2008a b) provide a more detailed 
discussion. The horizontal momentum, vertical momentum, 
continuity, and thermodynamic energy equations are 



dt 



-V$-/kx v + Pv 

9$ _ 1 
dp p 

V-V+— = 

op 

— = ± + — + Vt 

at Cp pcp 



(2.1) 
(2.2) 
(2.3) 
(2.4) 



where v is the horizontal velocity on constant-pressure sur- 
faces, Lo = dp/dt is the vertical velocity in pressure coor- 
dinates, $ is the gravitational potential on constant-pressure 
surfaces, / = 217 sin is the Coriolis parameter, Vl is the 
planetary rotation rate (In over the rotation period), k is the 
local vertical unit vector, q is the thermodynamic heating 
rate (Wkg"'), and T, p, and Cp are the temperature, den- 
sity, and specific heat at constant pressure. V is the hori- 
zontal gradient evaluated on constant-pressure surfaces, and 
d/dt = d/dt + y ■ V + ujd/dp is the material derivative. Cur- 
vature terms are included in v ■ Vv. Equation 12.41 is actually 
solved in an alternate form. 



d0 
dt 



9 q 
T c„ 



+ Ve 



(2.5) 



where 6 = Tip/ po)'^ is the potential temperature (a measure 
of entropy), k is the ratio of gas constant to specific heat at 
constant pressure, and /?() is a reference pressure (here chosen 
to be 1 bar, but note that the dynamics are independent of the 
choice of po). The dependent variables v, uj, $, p, 6, and 
T are functions of longitude A, latitude (p, pressure p, and 
time t. In the above equations, the quantities V^, Vj, and Vg 
represent additional source/sink terms beyond those described 
explicitly in the equations (see below). 
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Fig. 2. — Cubed-sphere grid at a resolution of C16 (i.e., 16 X 16 ele- 
ments per cube face), roughly equivalent to a global resolution of 64 X 32 in 
longitude and latitude. 

The MITgcm has been widely used in the atmospheric and 
ocean-science communities and has b een suc cessfully bench- 
marked against standard test cases jHeld & Suarez 1994. 
The model, with simplified forcing, has also shown success 
in reproducin g the global circulations of giant planets in our 
Solar System Eian & Showmanll2008l) . 

The MITgcm solves the equations on a staggered Arakawa 
C grid ( Arakawa & Lamb 1977) using a finite-volume dis- 
cretization. Two coordinate systems are supported: the stan- 
dard longitude/latitude coordinate system and a curvilinear 
coordinate system called the "cubed sphere," shown in Fig.|2l 
In the longitude/latitude system, the east-west (zonal) grid 
spacing converges to zero near the poles. This convergence 
implies, via the Courant-Fredricks-Levy (CFL) constraint, 
that the timestep must approach zero to maintain numerical 
stability. This difficulty can be surmounted using a polar filter 
that smooths the east-west variability of the dynamical fields 
near the poles, which increases the effective grid size and al- 
lows use of a finite timestep. In contrast, the cubed-sphere 
grid lacks the singularities at the poles and thereby sidesteps 
this problem. The price is a non-orthogonal grid and the exis- 
tence of eight special "corner points" (three of which can be 
seen in Fig. |2]l. Unlike the situation at the poles in a longi- 
tude/latitude grid, however, the coordinate lines in the cubed- 
sphere grid remain well separated until very close to the cor- 
ner points; depending on resolution, the size of finite- volume 
elements abutting the corners is typically one-half to one-third 
of that away from the corners. The upshot is that one can typ- 
ically use a longer timestep with the cubed-sphere grid than a 
longitude/latitude grid at comparable resolution. The simula- 
tions presented here adopt the cubed-sphere grid. 

The top boundary condition is zero pressure and the bot- 
tom boundary condition is an impermeable surface, which we 
place far below the region of interest. Both boundaries are 
free slip in horizontal velocity. We explore horizontal reso- 
lutions of CI 6, C32, and C64 (implying that each of the six 
"cube faces" has a resolution of 16 x 16, 32 x 32, and 64 x 64 
finite-volume elements, respectively); these roughly translate 
to global resolutions of 64 x 32, 128 x 64, and 256 x 128 in 
longitude and latitude. For a model with A^^. vertical layers, the 
bottom Nl-1 layers have interfaces that are evenly spaced in 



log pressure between ptop and the basal pressure; the top layer 
extends from a pressure of zero to ptop- In most simulations, 
we place the bottom boundary at 200 bars, safely below the 
active meteorology near the photosphere. We generally place 
Ptop at 0.2 mbar or 2 /xbar, using 40 or 53 layers, respectively; 
in both cases, this resolves each pressure scale height with 2.9 
layers. 

We a dopt the third -order Adams-Bashforth timestepping 
scheme (lDurraniri991h . which has favorable properties over 
the more commonly used leapfrog scheme. To maintain nu- 
merical stability, we apply a fourth-order Shapiro filter to the 
time derivatives of v and 6 at each timestep. This source/sink, 
which is represented by the terms 2?v and Vg in the govern- 
ing equations, is analogous to a hyperviscosity: it horizontally 
smooths grid-scale variations while leaving long-wavelength 
structures relatively unaffected. 

CFL constraints limit us to a short dynamical timestep; for 
our CI 6, C32, and C64 simulations, we generally use 50 sec, 
25 or 15 sec, and 6 sec, respectively. 

In our nominal simulations, we assume the planets rotate 
synchronously with their orbital periods, implying rotation 
periods of 2.2 days = 3.278 x lO-^s"') for HD 189733b 
and 3.5 days (n = 2.078 x 10"^ s"') for HD 209458b. Nev- 
ertheless, bec ause deviations from synch ronous rotation are 
possible (e.g. IShowman & Guillotll2002l) . we also ran some 
simulations with differing rotation periods. We adopt grav- 
ity and planetary radius of 21.4ms"^ and 8.2396 x lO^m 
for HD 189733b and 9.36ms-2 and 9.437 x lO^m for HD 
209458b. The equation of state is ideal gas. We fix Cp = 
1.3 X 10'*Jkg"'K"' and use k = 2/7, appropriate to apredom- 
inantly hydrogen atmosphere, for both planets. In most cases 
the initial condition contains no winds and adopts an initial 
temperature profile from a ID planetwide-average radiative- 
equilibrium calculation; 

As discussed by Showman et al.l ( l2008bl) and iGoodmaiil 
(l2009h . the mean winds that develop in an atmosphere de- 
pend in principle not only on forcing but also on damping 
processes that remove kinetic and/or potential energy. For 
solar-system planets, kinetic-energy loss occurs via turbu- 
lence, wave breaking, and friction against the surface (if any). 
In the deep interior of gas giants (p > 10^ bars), Ohmic 
dissipation may also provide an i mportant source of drag 
dKirk & Stevensonll 19871; iLiu et all 12008 ). although it is un- 
clear whether the atmospheric circulation on hot Jupiters cou- 
ples to such deep regions. Rigorously representing such fric- 
tional effects in numerical models is difficult, however; the 
turbulence associated with small-scale shear instabilities and 
wave breaking often has length scales much smaller (by up 
to several orders of magnitude) than can easily be resolved in 
global 3D numerical models. In Earth GCMs, these dissipa- 
tive processes are generally parameterized by adding to the 
equations empirical damping terms (e.g., a vertical diffusion 
to represent kinetic-energy losses by subgrid-scale shear in- 
stabilities and waves). A difficulty is that such prescriptions, 
while physically motivated, are generally non-rigorous and 
the extent to which they can be extrapolated to hot Jupiters is 
unclear. Moreover, simulations of the cloud-level circulation 
on Jupiter, Saturn, Uranus, and Neptune that lack frictional 
drag parameterizations have nevertheless shown success in 
reproducing the qualitative features of the observed flow on 
these planets (e.g., jc ott & Polvani 2007, 2008; Showmai3 
l2007l:[Lian & Showmanl2008l l2009). Therefore, while recog- 
nizing the possible importance that frictional processes could 
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play for hot Jupiters, for the present study we do not include 
frictional damping terms in the momentum equation (Eg. 12. Il l 
other than the Shapiro filter, which causes a diffusive damp- 
ing of small-scale wind structures. Once our simulations have 
spun up, this implies that, except for energy losses due to 
the Shapiro filter, the globally integra ted rate of production 
of available potential energy (APE; see lPeixoto &"Oora[T99l 
chapter 14) and its net conversion to kinetic energy (KE) be- 
come modest, with the potential-to-kinetic energy conversion 
in some regions of the atmosphere being counteracted by the 
kinetic-to-potential energy conversion in other regions.^ In a 
future study, we will explore the influence that plausible fric- 
tional parameterizations exert on the wind speeds and circu- 
la tion geometry. 

iGoo dman (2009) suggested that the heating associated with 
frictional dissipation should be included in the thermody- 
namic energy equation. Our study, like previous published 
hot- Jupiter studies, does not include this effect; this is equiv- 
alent to assuming that any frictional heating is small com- 
pared to the radiative heating. This is typically a reasonable 
assumption for atmospheres: the globally averaged rate of 
kinetic-energy production (and hence the rate at which that 
kinetic energy can be dissipated into thermal energy) is con- 
trolled by the thermodynamic efficiency of the atmospheric 
heat engine, which is typically much less than 100%, at least 
for atmospheres in the Solar System. In this situation the fric- 
tional heating is much less than the radiative heating. On 
Earth, for example, the globally averaged rate of frictional 
dissipation by the large-scale circulation is 2Wm"^, which 
is only a few percent of the latitude-dependent net radiative 
heating/cooling rate (Peixoto & Oort 1992). This implies that 
the frictional heating is only a small perturbation to the total 
heating rate. Because giant planets lack surfaces (which is the 
primary source of friction on Earth) is it plausible that the fric- 
tional damping, and the heating it causes, represents an even 
smaller effect on gas giants than terrestrial planets. Note that, 
to date, frictional heating has not been included in cloud-level 
atmosphere simulations of Jupiter, Saturn, Uranus, and Nep- 
tune. For hot Jupiters, given the considerable uncertainty aris- 
ing from other factors (e.g., uncertainties in opacities, compo- 
sition, possible presence of clouds/hazes, and possible non- 
synchronous rotation, all of which influence the circulation), 
it seems reasonable to neglect this effect for the present. 

2.2. Radiative transfer calculation 

We coupled the MITgcm to the plane -parallel radiative- 
transfer code of Marley & McKay (|1999|) . which is a state- 
of-the-art model that has been extensively used in ID in- 
vestigations of Tit an, the giant planets, brown dwarfs , 
and exoplanets (e.g . iMcKav et al.l ll989: M arlev et af]ll996l: 
Burrows et al."1997'; 'Marlev et al.''l999, 2002"; 'Fortney et ajj 
2005, 2006b a; Fortnev & Marley 2007; Fortney et al. 200i 
Saumon & Marlevj|20 08). We here use the code in its two- 
stream formulation, allowing treatment of the upward and 
downward radiative fluxes versus wavelength and height 
throughout the 3D grid. Multiple scattering is properly ac- 

' For example, our simulations develop a broad eastwai'd equatorial jet; 
near IR photosphere levels, where this jet transports air from day-to-night, the 
air generally flows in the same direction as the horizontal pressure-gradient 
force (i.e. v • V$ < 0), so potential energy is converted to kinetic energy. 
However, where the jet transports air from night-to-day, such a jet tends to 
flows against the horizontal pressure-gradient force (v ■ V'l' > 0), so kinetic 
energy is converted to potential energy. In a global mean, there is a large 
degree of cancellation between these competing effects. 



counted for. 

2.2.1. Opacities 

The code treats op acities using the correlated-k method 
( lGoodv& YundlI989 '), which is ideal for use in GCMs be- 
cause it is fast and accurate. In this approach opacities are first 
computed exactly at millions of individual frequency points 
at over 700 pressure-temperature (p-T) points, accounting 
for the influence of all molecular and atomic lines in our 
opacity database (iFreedman et alJ l2008l) . Line broadening 
is computed for each line of each species at each pressure- 
temperature point. The summed opacities at each p-T point 
thus account for both the relative abundances of each absorb- 
ing species as well as the line shapes appropriate for the par- 
ticular physical conditions. The spectrum of opacities is then 
divided into 30 frequency bins*^ as listed in Table [T] 

The opacities at each of the several tens of thousands of 
frequency point within each frequency bin (Table [T]) are then 
sorted by strength and described statistically thus enabling 
gaussian integration to be used to solve for the radiative fluxes 
through the model eight times, once for each probability inter- 
val within each spectral bin. The model flux within each of the 
30 frequency bins is the sum of the flux calculations weighted 
by the relative contribution of each statistical opacity inter- 
val. By using two intervals (one for the opacities weaker 
than 95% of the highest value and one for the strongest 5%), 
each with 4 gauss points (for 8 total gauss points), we resolve 
both the important strongest lines as well as the background 
continuum within each bin. This k-coefficient approach has 
been widely used in planetary atmospheres calculations for 
over two decades and is frequently employed in GCM calcu- 
lations. The method is subject to significant errors only when 
the distribution of opacity within a single spectral bin varies 
substantially with height in the atmosphere as might happen 
if one absorber is replaced by another with a very different 
opacity shape (for example a strong red slope within a bin re- 
placed by a strong blue opacity slope within that same bin). 
In practice such concerns are mitigated by careful choice of 
the opacity intervals. To compute the opacity at arbitrary p-T 
points we interpolate within our grid of k-coefficients. 

We chose the wavelength boundaries of our bins in an at- 
tempt optimize a number of competing needs. We wish to 
have highest spectral resolution (narrowest bins) in the opti- 
cal at the peak of the incident flux (to best resolve absorp- 
tion of incident energy) and in the near-infrared to similarly 
resolve the peak of the emergent flux. In addition, distinct 
bins should resolve important opacity peaks and valleys so 
that there is a relatively flat opacity distribution within each 
bin. This minimizes the dynamic range that must be covered 
by the 8 k-coefficient Gauss points within each bin. Exam- 
ples include the water opacity windows in the near-infrared, 
the CO bandhead at 2.3 fim, and TiO in the optical. We also 
wish to minimize the total number of bins while still covering 
the entire spectral range from the UV to the thermal infrared. 
Our choice of bins is listed in Table[T]and depicted graphically 
in Fig.|3] 

The opacities themselv es are computed from a compre- 
hensive opacity database ( IFreedman et al.l 12008,) , calculated 
assuming local thermodynamic and chemical equilibrium 

^ In recent ID applicatio n s of ou r code we typically use 196 intervals (e.g. 
IFortnev et al.l2OO5I. E006bl l3. l2008ll . Extensive testing has revealed that, for a 
specified T(p) profile, the net radiative flux calculated using 30 bins differs 
less than 1% from that calculated using 196 and is substantially faster in the 
GCM enviroiunent. 
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Fig. 3. — Opacity (in cm^ molecule"' ) versus wavelength at p = 1 bar and 
T = 1 000 K (red) and 2000 K (blue) assuming solar composition (modified for 
rainout where appropriate) including TiO and VO. In the mid-IR, much of the 
structure results from H2O, CO, and (at 1000 K) CH4; the two prominent Na 
lines are evident in the visible in the 1000-K case. Significant additional 
visible-wavelength opacity occurs in the 2000-K case due to absorption by 
TiO and VO. Wavelength boundaries for the 30 opacity bins in our GCM 
runs are shown with vertical orange lines and listed numerically in Table[T] 



TABLE 1 
Opacity Bins 



wavelength (^m) wavelength (/^m) 



324.68 


46.00 


46.00 


20.00 


20.00 


10.40 


10.40 


6.452 


6.452 


5.220 


5.220 


4.400 


4.400 


3.800 


3.800 


3.288 


3.288 


2.989 


2.989 


2.505 


2.505 


2.170 


2.170 


2.020 


2.020 


1.777 


1.777 


1.593 


1.593 


1.497 


1.497 


1.330 


1.330 


1.197 


1.197 


1.100 


1.100 


1.005 


1.005 


0.960 


0.960 


0.910 


0.910 


0.860 


0.860 


0.785 


0.785 


0.745 


0.745 


0.675 


0.675 


0.612 


0.612 


0.572 


0.572 


0.495 


0.495 


0.400 


0.400 


0.261 



Note. — Opacity bins used in the 
radiative-transfer calculation. Each row 
lists the bounding wavelengths for one 
bin. 

jLodders & FeglevI l2002l l2006h as a function of p, T, and 
wavelength. We perform simulations with opacities corre- 
sponding to 1, 5, and 10 times solar abundances (appropri- 
ately modified for rainout of condensed species whenever 
condensation sho uld occur). For HP 189733b, the infer red 
water abundance (iTinetti et al.ll2007t ISwain et all l2008bl) is 
consistent with near-solar metallicity for both wa ter an d car- 
bon dShowman 2008). For HD 209458b, Barman ( 2007h 's fits 
to transit data suggest abundances not exceeding a few times 
solai", while Sing et al.. (2008) 's fits to transit data suggest two 
times solar sodium abundance. 



TiO and VO deserve special mention. For our sim- 
ulations of HD 189733b, we use opacity databases that 
exclude TiO and VO, consistent with secondary-eclipse 
data dCh arbonneau et al. 2008), ID radiative-transfer studies 
(iBarmanI 12008: Knut son et ah 2009), and evolution models 
suggesting that TiO and VO should sequester deeper than 
~ 10- 100 bars where solid clou ds incorporating Ti and V 
form (e.g. iFortnev et all l2006tl) . For HD 209458b how- 
ever, the detection o f a hot stratosphere (iKnutson et al.ll2008l : 
iBurrows et al. 2007), possibly associated with TiO and VO 
dHubeny et al. 2003: Fortnev et al . 2008) suggest that gaseous 
TiO and VO may exist in this planet's atmosphere. For HD 
209458b, therefore, we run some cases that include gaseous 
TiO and VO and other cases that exclude them. 

For simplicity, we neglect cloud opacities. This appears 
to be reasonable for HD 189733b, at least for emitted ther- 
mal radiation, since absorptio n features observed in transit 
spectra long ward of ^ 1 .5 /zm (iTinetti et alj2007l : ISwain et alj 
l2008b: Sing et al.]l2008h indicate that gas rather than particles 
dominates the opacity along the slant transit pat h. The sit- 
uation is less clear at visible wavelengths ( Pont et al.l 120081 : 
iRedfield et al.ll200^ . but submicron-sized haze particles are a 
possibility. In any event, because the slant optical d epth ex- 
ceeds the normal optical depth by a factor of ~ 50 dFortnevI 
l2005h . the data are nevertheless consistent with normal cloud 
optical depths <C 1. Fewer constraints exist for HD 209458b. 
Some authors hav e suggested silicate ab sorption features in 
emission spectra ([Richardson et al J [20070 . but debate exists 
(ISwain et alJl2008alK 

2.2.2. Flux Calculation 

We solve for the net flux through each level of a specified 
local atmospheric structure using two separate techniques: 
one for the incident stellar and one for the emitted planetary 
thermal radiation. This subdivision is common in planetary- 
atmosphere studies since the incident radiation is dominated 
by the direct stellar beam accompanied by a weaker but more 
uniform scattered flux while the thermal emitted radiation is 
comparatively much more uniform over a hemisphere, thus 
suggesting separate flux solution methods. For solar-system 
atmospheres the solar and planetary fluxes usually occupy dis- 
tinct spectral intervals. However for hot extrasolar planets the 
incident and emitted fluxes substantially overlap near 1 /im 
and thus the net flux within a layer is the sum of the net stel- 
lar and thermal fluxes, which are computed by differing tech- 
niques within the same spectral interval. 

For the incident stellar flux w e follow iMcKav et alj 
(119891) and iMarley & Mc"Kay[ ([1999 ) and emp loy the delta- 
discre te ordinates method (see Appendix I of iMcKav et alj 
(Il989l) ). We adopt stellar spectra from Kurucz (see 
[http : //hurucz .harvard . edu/s tars . html ) and assume that HD 
189733b and HD 209458b have distances of 0.0313 and 0.046 
AU from their stars, respectively. Obliquities and orbital ec- 
centricities are assumed zero. Stellar fluxes are computed 
only in those bins lying shortwards of 6 /im. 

For the thermal flux we follow the approach used in our 
brown dwarf work and employ the " two-stream source func- 
tion" technique of iToon et al.l (Il989i) which was specifically 
developed for rapid, accurate computation of radiative heat- 
ing rates in inhomogeneous atmospheres. This approach 
is well suited for atmospheres with scattering cloud layers, 
which we will eventually include. The present implementa- 
tion does include Rayleigh scattering, although this is gener- 
ally unimportant for most of the thermal bins. The accuracy 
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of this technique for vari ous hmiting cases is discussed in de- 
tail by iToon et alJ (Il989h . We model thermal radiation from 
the planet from 0.26-325 iim. Extensive experience with this 
technique employed with k-coefficient opacities in the context 
of brown dwarf atmospheres has shown that integrated fluxes 
usually agree to much better than 0.1% with fluxes computed 
using detailed line-by-line calculations. 

In most ID radiative-transfer calculations, an iteration is 
performed to drive the solution into radiative equilibrium. 
Here, the 3D solution is not in radiative equilibrium, and no 
iteration is needed. Instead, at each timestep, we pass a driver 
subroutine the most recent 3D temperature field, consisting 
of numerous individual T(p) columns, one per element in the 
horizontal grid. We loop over all these columns and, for each, 
we call the radiative-transfer solver once (without iteration) to 
determine the upward and downward fluxes versus frequency 
and pressure throughout that column. For each column, we in- 
clude the appropriate incident starlight; columns on the night- 
side emit radiation but receive no starlight, while columns on 
the dayside receive starlight along an angle 9 from the local 
vertical. The cosine of this angle, /i, varies from 1 at the sub- 
stellar point to at the terminator Columns closer to the sub- 
stellar point thus receive greater stellar flux, and this flux pen- 
etrates to greater pressures because of the shorter atmospheric 
path lengths. 

2.3. Coupling dynamics and radiative-transfer 

Once we calculate the wavelength-dependent radiative 
fluxes, we sum them to obtain the net vertical flux everywhere 
over the 3D grid at that timestep. We then obtain the thermo- 
dynamic heating rate q (Eq. 12.41 ) as follows. If an atmospheric 
layer absorbs more (less) radiation than it emits, the layer ex- 
periences heating (cooling). The thermodynamic heating rate 
per volume ( Wm"-') is thus simply -dF/dz, where z is height 
and F is the net, wavelength-integrated flux. Expressed in 
pressure coordinates using hydrostatic balance, the heating 
rate per mass, q, becomes 

dF 

q = 8jr (2.6) 
op 

This is then used to drive the dynamics via Eq. l2.5l (see Fig.[T]i. 
Note that this expression neglects the heating contribution due 
to the horizontal divergence of the horizontal radiative flux, 
which is at least a factor of ^ 30 less than that in Eq. |2.6| for 
conditions relevant to hot Jupiters. 

We use the same vertical gridding for the dynamics and 
radiative-transfer calculations. In the SPARC/MITgcm, the 
vertical grid is staggered such that T, 9, and horizontal wind 
are defined within the layers, while certain other quantities 
(such as vertical velocity) are defined at the interfaces be- 
tween layers. Because q is used to update 9 (Eq. |2.5t , q should 
be computed within the layers so that it has the same vertical 
positioning as 9. To calculate q we thus evaluate dF /dp by 
finite differencing the fluxes and pressures between the inter- 
faces that over- and underlie a given layer. This maximizes 
accuracy while positioning dF /dp within the layer Note that 
we use the same layer spacing in the dynamics and radiative 
transfer calculations; ID radiative transfer tests performed of- 
fline suggest that this vertical gridding is sufficient to resolve 
the radiative transfer 

As a test, we ran a suite of simulations where we turned off 
the dynamics. In these cases, the temperature profiles relaxed 
toward the expected radiative equilibrium (as defined by solu- 
tions generated with the radiative-equilibrium ID version of 



the radiation code) for the specified value of p. This demon- 
strates that the calculation of the thermodynamic heating rate 
(Eq. 12.61 ) functions properly. 

In most SPARC/MITgcm simulations, we update the radia- 
tive fluxes less frequently than the dynamical timestep. This 
is standard practice in GCMs and helps maintain computa- 
tional efficiency. For our simulations without TiO and VO, 
we found that a radiative step of 500-lOOOsec was generally 
sufficient. For simulations with TiO and VO, however, a much 
shorter radiative timestep of 100 sec or less was necessary to 
maintain numerical stability. For computational expediency, 
in some simulations, particularly our high-resolution cases, 
we perform radiative transfer only on every other grid column 
and interpolate the fluxes and heating rates in between. 

SPARC/MITgcm output includes the emergent fluxes ver- 
sus wavelength and position everywhere over the globe. How- 
ever, because we run the GCM with only modest spectral res- 
olution, we generally recalculate disk-averaged light curves 
and spectra offline at higher spectral resolution using the 
3D temperature fields following the procedure described in 
Fortnev et al.. (.2006a) and Showman et al. (2008a). This pro- 
cedure ensures that limb darkening is properly treated. 

3. RESULTS: SYNCHRONOUSLY ROTATING HD 
189733B 

3 . L Circulation regime 

In our simulations, the imposed day-night irradiation gradi- 
ent rapidly leads to the development of large day-night tem- 
perature contrasts and winds reaching several km sec"' . A sta- 
tistically steady flow pattern is achieved at pressures < 1 bar 
after less than 100 Earth days of integration. Figures |4]-(5] 
illustrate the behavior for our nominal, synchronously rotat- 
ing simulation of HD 189733b using solar opacities without 
TiO and VO. Figure |4] shows the temperature (colorscale) and 
winds (arrows) over the globe on three different pressure lev- 
els, while Fig.|5]depicts the zonally averaged zonal wind^ ver- 
sus latitude and pressure. 

As in our previous simulations with simplified forcing 
(Showman & Guillot)|2002t iCooper & Showmanll200l l200a 
Showman et al. 2008a|), the flow becomes dominated by a su- 
perrotating (eastward) equatorial jet extending from latitudes 
of approximately 30°N to 30°S. Vertically, the jet remains 
coherent from the top of the model (0.2 mbar) to the ^ 10- 
bar level and reaches peak zonal-mean zonal wind speeds of 
3.5kmsec"' at pressures of 10-100 mbar (Fig. |5]l. The lat- 
itudinal width and longitudinal structure of the jet depend 
strongly on altitude. Longitudinal variability in the jet speed 
is small at pressures exceeding ^ 200 mbar but increases with 
altitude and exceeds 2 km sec"' at the top of the model, with 
the fastest eastward winds occurring in the hemisphere west 
of the substellar point and the slowest winds occurring near 
and east of the substellar point. The transition between these 
regions narrows with altitude, and at pressures less than a few 
mbar, the equatorial flow appears to pass through a hydraulic 
jump as it approaches the substellar point from the west (e.g., 
Fig.|4] top panel), leading to a near-discontinuity in the zonal 
and meridional wind speed at longitudes of -50° to -80° (de- 
pending on latitude). Accompanying this structure are strong 
downward velocities, and the adiabatic compression associ- 

' In atmospheric dynamics, tlie terms zonal and meridional wind refer to 
the east-west and north-south wind, respectively (with eastward and north- 
ward being positive); a zonal average refers to an average of any quantity in 
longitude. 
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FiG- 4. — Temperature (colorscale, in K) and winds (an'ows) for nomi- 
nal HD 189733b simulation with solar abundances and no TiOA'O. Panels 
show flow at 1 mbar (top); 30 mbar, corresponding to an approximate pho- 
tosphere level in the mid-IR (middle); and 1 bar (bottom). Resolution is C32 
(roughly equivalent to a global horizontal resolution of 128 X 64 in longitude 
and latitude) with 40 vertical layers- Substellar point is at longitude, latitude 
(0° , 0° )- Dayside is region between longitudes -90° and 90° . Nightside is at 
longitudes -180° to -90° and 90° to 180°- 

ated with this descent leads to a warm chevron-shaped struc- 
ture in the temperature (Fig.|4] top panel). 

At higher latitudes (poleward of 45°), the zonal-mean 
zonal wind is westward at pressures < Ibar (Fig. |5]l. Nev- 
ertheless, this high-latitude flow exhibits complex structure, 
with westward flow occurring west of the substellar point, 
eastward flow east of the substellar point, and substantial 
day-to-night flow over the poles. These polar flows indicate 
the importance of performing global simulations and would 
not be properly captured if the poles were removed from 
the computation dornain, a s done by some previous authors 
dPobbs-Dixon & Linll2008l) . 

On average, the dayside is hotter than the nightside, but 
the dynamics distorts the temperature pattern in a complex 
manner (Fig. [Hi. Perhaps most importantly, the hottest re- 
gions do not occur at the substellar longitude; instead, ad- 
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Fig . 5 - — Zonal-mean zonal winds for our nominal HD 1 89733b simulation 
at solar abundance. This is the same simulation as in Fig.|4] Scalebar gives 
speeds in m s"' . 



vection associated with the equatorial jet shifts the hottest 
regions downwind (eastward) of the substellar point by 
~ 50°. This point has been previously emphasiz ed by 
Showman & Guillot (2002'), 'Coope r & Showrnanl 12005), and 
Showman etal. (2008a). Moreover, the coldest regions do 
not occur at the equator but instead within two broad gyres 
centered at latitudes of ±40-50° and longitudes ~ 60-80° 
east of the antistellar point — a phenomenon not seen in our 
previous simulations with simplified forcing (IShowman et al.l 
2008a). The horizontal wind speed is almost zero near the 
center of these gyres, so air parcels trapped there have long 
residence times on the nightside. This allows them to expe- 
rience a large temperature drop due to radiative cooling. In 
contrast, air within the equatorial jet has only a short resi- 
dence time (typically ^ 1 day) on the nightside because of the 
fast jet speeds, leading to only modest temperature decreases 
via radiative cooling. Interestingly, the temperature structure 
shows substantial longitudinal variability even at pressures as 
great as 1 bar. 

As co mpared to our previous simulations with simplified 
forcing (Showman et all l2008al) . our current HD 189733b 
simulations exhibit modest lateral temperature contrasts. The 
horizontal temperature differences reach ^ 450 K at 1 bar and 
^ 750 K at 1 mbar (Fig.|4|i. In contrast, in our previous simu- 
lations forced by Newtonian heating/cooling (Showman et ai] 
|2008a), the day-night temperature differences reached nearly 
900 K at 100 mbar and 1 000 K at 10 mbar. The smaller values 
in our present simulations, which can be attributed to our us- 
age of realistic radiative transfer, have major implications for 
light curves and spectra (§3.2). 

Notice that the global-scale temperature structure exhibits 
significant vertical coherency throughout the observable at- 
mosphere (Fig. |4]i. Although the detailed structure varies 
between levels, the hottest regions lie east of the substellar 
point throughout, with the longitudinal offset of the hottest 
region varying only modestly between pressures of 1 bar 
and 1 mbar. Likewise, the locations (though not the shape) 
of the coldest regions also maintain coherency across this 
pressure range (blue regions in Fig. |4|l. At first glance, 
this vertical coherency is surprising, because idealized ra- 
diative calculations have shown that the radiative time con- 
stant s hould vary by o r ders of magnitude over this pressure 
range (liro et all 120051: iFortnev et all 120081: IShowman et"!!] 
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l2008ab . At pressures where the radiative time constant is 
comparable to the time for wind to advect across a hemi- 
sphere, one expects a signific ant offset o f the hottest regions 
from the substellar point ( Sho wman & G uillot 2002). How- 
ever, at low pressures where the expected radiative time con- 
stants are much shorter than plausible advection times, one 
expects the temperature patterns to track the stellar heat- 
ing, w ith the hottest regio n occurr ing close to the sub stellar 
point (ICooper & Showrna n 2005; Knutson et alJl200"7l) . In- 
deed, precisely this height-dependent pattern is seen in pub- 
lished 3D simulations that force the flow with a simplified 
Newtonian heating/cooling scheme, which relaxes the tem- 
perature toward the radiative-equilibrium temperature profile 
over the expected radiative timescale (Showman et al. 20082 
iFortnev et a l. '2006a: I Cooper & Showmanll2005L |2006). 

So what causes the vertical coherency in our current simu- 
lations? The simple arguments described above — in which 
the temperature approaches radiative equilibrium if radiation 
times are less than advection times — implicitly assume that 
the radiative-equilibrium temperature profile can be indepen- 
dently defined and that it has a structure reflecting that of the 
insolation, with the greatest radiative-equilibrium temperature 
at the substellar point. However, this argument neglects the 
fact that, in real radiative transfer, the radiative-equilibrium 
temperature profile /fie// depends on the dynamical response 
and can involve radiative interactions between different lev- 
els.'" To illustrate, suppose the small heating rates at ~ Ibar 
lead to a hot region shifted east of the substellar point (as 
seen in Fig. HI bottom panel). Upwelling infrared radiation 
from these hot, deep regions will warm the entire overlying 
atmosphere, leading to a temperature pattern at low pressure 
that has similar spatial structure as that at higher pressure. 
These effects were ignored in our previous s tudies adopt- 
ing Newtonian cooling ([Sh owman et al. 2008aL IFortnev et aP 
I2006al: ICooper & Showmanii2005>. i2006 ). but they are self- 
consistently included here and can explain the vertical co- 
herency in Fig. m despite the large expected vertical varia- 
tions in radiative time constant. Nevertheless, our 5 and 10 
times solar HD 189733b simulations (and especially the HD 
209458b simulations with TiO and VO opacity to be discussed 
in §5) exhibit less vertical coherency than shown in Fig.|4] 

Despite these quantitative differences, our current sim- 
ulations lie within the same basic dynamical regime as 
our previous 3D simulations driven by Newtonian coo ling 
dShowman et al.ll2008al: ICooper & Showmanll2005l |2006'). In 
all these cases, the flow exhibits a broad, eastward equatorial 
jet with westward flow at high latitudes; eastward displace- 
ments of the hottest regions from the substellar point (at least 
over some range of pressures); and a gradual transition from 
a banded flow at depth to a less-banded flow aloft. More- 
over, in all these cases, the flow structures have horizontal 
lengthscales comparable to the planetary radius, consistent 
with the large Rhines scale and Rossby deformation radius for 
fliese planets (IShowman & Guillotil2002t iMenou et al.ll2003t 
IShowman etalj|2008bl) . 

Figure |6] top panel, illustrates the diversity of vertical tem- 
perature profiles that occur Red (blue) profiles lie equator- 
ward (poleward) of 30° latitude. This is for our nominal HD 
189733b simulation with solar abundance; qualitatively sim- 

A more rigorous way of stating this is that one cannot define a radiative 
equilibrium temperature structure in isolation from the dynamics. In that 
case, the comparison-of-timescales argument fails and one must solve the 
full radiative-dynamical problem, as we are doing here. 
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Fig. 6. — A selection of profiles of temperature (top) and potential temper- 
ature (bottom) versus pressure in our nominal, solar-abundance HD 189733b 
simulation (the same simulation as in Fig.|4). Red (blue) profiles are equator- 
ward (poleward) of 30° latitude. 



ilar patterns occur for 5 and 10 times solar. Key points are 
as follows. First, as expected, equatorial regions are on av- 
erage warmer than the high latitudes at pressures less than 
a few bars. At low pressures, longitudinal variation is com- 
parable to the latitudinal variation. Second, the temperature 
declines smoothly with altitude from ^ 3 bars to ^ lOmbar 
This has important implications for spectra and light curves, 
which originate within this layer. Third, dynamics has modi- 
fied the deep stable radiative layer from 10-100 bars, leading 
to significant latitude variation of both temperature and static 
stability, with warm poles and a cold equator. 

Our HD 189733b simulations remain convectively stable 
everywhere throughout the atmosphere. This is illustrated in 
Fig. |6] bottom panel, which shows the potential temperature 
6 versus pressure for the same profiles depicted in the top 
panel. Potential temperature is a measure of entropy; an at- 
mosphere that is neutrally stable to convection has 9 constant 
with height, whereas 9 increases wi th altitude in a statically 
stable atmosphere (e.g. 'S albvlll99d pp. 166-172). As can 
be seen in Fig. |6] bottom panel, 9 increases with altitude in 
all the profiles. This is true even on the nightside. Nightside 
cooling (which tends to be stronger at low pressure than high 
pressure) reduces the static stability of the nightside profiles 
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Fig. 7. — Emergent flux density (ergs s"' cm^^Hz"') from our nominal, 
solar-abundance simulation of HD 189733b at six orbital phases. Black, 
nightside, as seen during transit; red, 60° after transit; green, 120° after tran- 
sit; dark blue, dayside, as seen during secondary eclipse; light blue, 60° after 
secondary eclipse; and magenta, 120° after secondary eclipse. The key in 
the top right corner is color-coded with the spectra to illustrate the sequence. 
Thin dotted black lines at the bottom of the figure show normalized Spitzer 
bandpasses and the letters at the top show locations of the H, K. L, and M 
bands. This is the same simulation as in Fig.O For comparison, the dot- 
ted curve is a spectrum from a ID planetwide average radiative-equiUbrium 
model. 

relative to the dayside profiles, but in our simulations this ef- 
fect is insufficient to generate a detached convection layer on 
the nightside. 

Simulations of HD 189733b performed using opacities cor- 
responding to 5 and 10 times solar abundances exhibit simi- 
lar behavior to our solar case depicted in Figs. |4}{6] The pri- 
mary difference is that, compared to the solar case, the 5 and 
10-times-solar cases have slightly warmer daysides, slightly 
cooler nightsides, and smaller eastward offsets of the hot re- 
gion from the substellar point, especially at low pressure. At 1 
bar, in all three cases, the approximate centroid of the hottest 
regions lie ~ 80° east of the substellar point. At 1 mbar, how- 
ever, the offset drops to ^ 50° in the solar case but only ~ 20° 
in the 5 and 10-times-solar cases. These differences can be 
attributed to the greater opacities in the 5 and 10-times-solar 
cases, which lead to greater heating rates and hence shorter 
radiative time constants. Despite the differences, we empha- 
size that all three cases exhibit extremely similar qualitative 
temperature and wind patterns. 

3.2. Spectra and light curves 

We now turn to spectra and light curves. Because our sim- 
ulations couple dynamics to radiative transfer, our model out- 
put includes the emergent spectral flux versus wavelength 
everywhere on the model grid at each timestep. How- 
ever, because we use a relatively small number {N\ = 30) of 
correlated-k wavelength bins, any spectra calculated from this 
output would have coarse spectral resolution. Thus, we in- 
stead re-calculated the emergent fluxes offline at greater spec- 
tral resolution (Nx = 196). To do so, we took the 3D temper- 
ature output at a given time and ran each vertical T{p) col- 
umn through our ID spectral solver to calculate upward and 
downward fluxes versus wavelength, taking care to use the 
same opacity database and stellar model as used in the orig- 
inal 3D simulation. The resulting emergent fluxes are essen- 
tially identical to those generated by the 3D model except that 



they have higher spectral resolution. We then calculate disk- 
integrated S2ectra_andjiehtcurves us ing the procedures de - 
scribed in lShowman et al.l (l2008ah and lFortnev et al.l (l2006al) . 
Planetary limb darkening (or brightening), as viewed by the 
observer, is accounted for 

Figure |7] shows the resulting spectra for our solar- 
abundance HD 189733b simulation at six orbital phases. Be- 
cause the simulated temperatures decrease with altitude in the 
relevant pressure range (s^ 0.01-1 bar), spectral features are 
seen in absorption (rather than emission) throughout the or- 
bit. Near the time of transit {black curve), the nightside faces 
Earth, with its cool temperatures and large vertical tempera- 
ture gradients (Fig.|6]l. This leads to relatively low fluxes with 
deep absorption bands of H2O and CH4. The hotter dayside 
(solid dark blue curve) has a smaller temperature gradient, 
leading to greater fluxes and shallower absorption bands. In- 
terestingly, the fluxes 60° of phase before transit (magenta) 
are smaller than those occurring at transit, while the fluxes 
60° of phase before secondary eclipse (green) exceed those 
occurring immediately before/after the secondary eclipse it- 
self. This phase offset results directly from the fact that the 
hottest and coldest regions are shifted eastward from the sub- 
stellar and antistellar points, respectively (see Fig.|4]i. 

Figure [8] shows light curves calculated in Spitzer band- 
passes for our HD 189733b simulations with opacities cor- 
responding to solar (top panel) and five times solar (bottom 
panel) abundances, respectively. Black, red, green, dark blue, 
light blue, and magenta show the simulated light curves at 
3.6, 4 5, 5.8, 8, 16, and 24 ^m. Overlaid are the lKnutson et all 
(TWf) 8 -/im light curve in small blue dots and the binned 
Knutson et al. (2009) 24-/im light curve in small magenta di- 
amonds, along with the Sp itzer secondary-echpse d epths from 
ICharbonneau et al.l (l2008h and iDeming et aH (l2006i) in large 
diamonds. 

Overall, our simulated light curves (Fig. |8]l compare favor- 
ably to the observed ones. We are able to reproduce the mod- 
est day-night flux variation seen in the observations at both 8 
and 24 /im; this contrasts wi th our earlier simulatio ns using 
Newtonian heating/cooling dShowman et al.l l2008al) . which 
greatly overpredicted the day-night flux variation. In our 
cuiTent solar-opacity simulations, the ratio of maximum-to- 
minimum flux (within a given Spitzer channel) ranges from 
1.4 to 3.5 depending on wavelength, while in the five-times- 
solar case it ranges from 1.6 to 4.1, with the greatest flux ra- 
tios occuiTing at 3.6 /im and the smallest at 16 and 24 /im. 
Likewise, our simulated light curves reach their peak flux be- 
fore the secondary ec lipse, a feature shared by bo th the 8- and 
24-/im Ught curves (iKnutson etal.ll2007l l2009l) . In the so- 
lar case (top panel), the offsets are close to 50°, whereas at 
five times solar (bottom panel), the offsets range from ~ 26° 
at 24 /im to ^ 42° at 3.6 /im. In the simulations, this phase 
offset results directly from the eastward displacement of the 
hottest and coldest regions from the substellar and antistel- 
lar points, respectively (Fig. |4|i. This phenomenon also oc- 
curred in our previous simulations for ced with Newtonian 
heating/cooling (Showman et al. 2008a; Cooper & ShowmanI 
'2005; Showman & Guillot 2002). 

We emphasize that the simulated light curves in Fig. [8] are 
not fits to the observations; beyond choosing the metallicity, 
no tuning of any kind was performed. Instead, Fig. [8] dis- 
plays the natural interaction of radiation and dynamics as re- 
solved by the model. Indeed, by explicitly representing both 
the dynamics and the radiation, our goal here is to eliminate 
the tunable knobs that have been used to parameterize dynam- 
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Fig. 8. — Light curves versus orbital pliase calculated in Spitzer band- 
passesfor HD 189733b. Top and bottom panels show light curves for our 
simulations using solar and five times solar abundances in the opacities, re- 
spectively. Within each panel, moving from bottom to top, the light curves 
are for wavelengths 3.6 fim (black), 4.5 /im (red), 5.8 fim (green), 8 ^m (dark 
blue), 16 /^m (light blue), and 24fim (magenta), respectively. Overplotted is 
the Spitzer Sfim light curve from Knutson et al. (2007) in dark blue points 
and the binned 24pm light curve from Knutson et al. (2009) in small ma- 
genta diamon ds. Spitzer secondary-eclipse depths from Charbonneau et al. 
(2008) and De ming et alj 12006) are plotted at 180° phase in large diamonds, 
with wavelengths color-coded as described above. Both simulations have a 
resolution of C32 with 40 layers (top panel is same simulation as in Fig.|4). 

ics and/or radiation in some previous studies. 

Nevertheless, there exist some important discrepancies be- 
tween the simulated and observed light curves. First, we do 
not reproduce the flux minimum that occurs ^ 50° after transit 
in the observed 8-/xm light curve. If real, this feature suggests 
the exi stence of a local cold region to the west of the antistellar 
point (iKnutson et al.ll2007t) . However, this flux minimum i s 
absent in the Spitzer 24- /im light curve (IKnutson et al.ll2009l) . 
and analysis of these data suggest instead that the minimum 
flux region actually lies east of the antistellar point, which 
would be qualitatively consistent with our simulations. This 
hemisphere of the planet's surface is not well-resolved by the 
existing light curves, which cover only half an orbit; more 
data are neede d to resolve the issue of where the true flux 
minima lie (Kn utson et al.ll2069h . 

Second, the phase offsets in our simulations are somewhat 
too large, especially in our solar-opacity case. The observed 
offsets of the flux peak are 16° ± 6° at 8 ^lm and 20-30° at 
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Fig. 9.— Con tribution functions (e.g. IChamberlain & Huntei] I1987t 
IKnutson et alj|2009l) calculated using our ID radiative transfer model for a 
generic cloud-free pL-class planet without atmospheric TiO and VO (left) 
and a generic pM-class planet with atmospheric TiO and VO (right). Both 
are for dayside conditions, and both assume solar metallicity with equilib- 
rium chemistry. Contribution functions are calculated for various Spitzer 
broadband filters (black short dashed, red solid, green dashed-dotted, blue 
dashed-triple-dotted, and pink long-dashed curves for 3.6, 4.5, 5.8, 8, and 
24 pm, respectively), K band (orange dotted curve), and the Kepler band at 
450-900 nm (black soUd curve). For clarity some of the curves have been 
normalized to 0.5 or 0.75 rather than 1. 

24 //m. In contrast, in the solar abundance simulation (top 
panel), the phase offsets are close to 50° at all Spitzer band- 
passes. The agreement is better in 5-times-solar case, where 
the simulated offsets are 30° at 8 /im and only 26° — per- 
fectly consistent with the observed offset — at 24 /im. On the 
other hand, at 24 /im, the solar case provides an overall better 
match to the magnitudes of the light curve flux values. 

We now turn to the secondary-eclipse observations (large 
diamonds in Fig. |8]l. We match reasonably well the eclipse 
depths at 4.5, 5.8, 8, 16, and 24 /im, especially with our five- 
times-solar-metallicity model, although it would appear that 
our solar model is slightly too cool. The greatest discrepancy 
in both models is that we underpredict the 3.6- fim secondary 
eclipse depth, by factors of 2.1 and 1.8 in our solar and five- 
times-solar cases, respectively. Our radiative-transfer calcu- 
lations suggest that 3.6 /im samples pressures of ^ 0.1-1 bar, 
deeper than any other Spitzer bandpass (see contribution func- 
tions in Fig.|9] left panel). One possibility is that our model is 
too opaque at this bandpass; a lower opacity would allow pho- 
tons to escape from greater pressures, where temperatures are 
hotter (Fig. |6]l. Fits to the spectra of T dwarfs with the same 
atmospheric chemistry and radiative transfer and similar ef- 
fective temperatures also suggest too much model opacity in 
this IRAC bandpass (Geb alle et al.ll2009i) . Alternatively, our 
model could simply be too cold in the 0.1-1-bar region, at 
least on the dayside. Interestingly, published ID models of 
HP 189 733b also have difficulty matching the 3.6 /im point 
dBarm an 2008; Knutson et al. 2009). A fuU Ught curve at 
3.6 /im, possible with warm Spitzer, would be invaluable in 
constraining this deep layer of the atmosphere. 

Figure [To] reiterates these points by displaying the planet- 
to-star flux ratio versus wavelength. Our solar and five-times- 
solar cases are depicted in red and green, respectively. Photo- 
metric data are shown in blue. The reasonable agreement at 
4.5, 5.8, 8, 16, and 24 /im is evident, as is the discrepancy at 
3.6 /im. Interestingly, this disrepancy is worse in our model 
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Fig. 10. — Planet-to-star flux ratios for our HD 189733b simulations at 
the time immediately before/after secondary eclipse. Red and green de- 
pict our solar and 5 t imes solar cas es, respectively. Triangle shows the 
16/im IRS point from 'Deming et al.' f2006). Blue squares show data from 
ICharbonneau et al. (2008), including their re-analysis of the Deming et al] 
(2006) data at 16/jm. Circles at 8 and 24 ^m give secondary-eclipse depths 
obtained from the light curves of Knutson et al. (2007, 2009). Line at 2.2 
(K band) gives the upper Umit from Barnes et al. (2007). Grey points give 
the Spitzer IRS spectrum from Grillmair et al. (2007). Inset: Planet-to-star 
flux ratios on the nightside for these same models. Thin blue circles show 
the nightside flux ratios at 8 and 24 /jm as obtained from the light curves of 
iKnutson et al. (2007, 2009). Thick blue circles are those same data corrected 
for the effect of starspots. 
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Fig. 1 1 . — Effective temperature versus orbital phase for the Earth-facing 
hemisphere of our nominal solar-opacity simulations of HD 189733b (thick 
curve) and HD 209458b (thin curve). 

than in published ID calculati ons. On the other hand, un- 
like c urrent ID dayside models (lBarmanll2008l: iKnutson et aTl 
20091). we ar e able to satisfy the upper hmit at 2.2 /xm from 
Barnes et alj (12007). On the nightside (see inset), we fare rea- 
sonably well against the 8-/im and 24-/im planet-to-star flux 
ratios from the light curves {blue circles). 

Figure [TT] shows the total luminosity versus orbital phase 
for our solar-opacity simulation of HD 189733b expressed as 
an effective temperature. This was calculated by integrating 
the planet's spectrum (Fig. |7J over all wavelengths to obtain 
a flux, dividing by the Stefan-Boltzmann constant, and taking 
the fourth root to obtain a temperature. The effective tem- 
perature reaches minima and maxima of ^ 1050 and 1250 K, 
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Fig. 12. — Sensitivity of light curves, calculated in Spitzer bandpasses, to 
model resolution (top) and integration time (bottom). Top panel shows HD 
189733b simulations at horizontal resolutions of C64 (solid), C32 (dashed), 
and C16 (dotted), all at 301 Earth days of integration time. Bottom panel 
shows HD 189733b C16 simulations at 104 (dotted), 995 (dashed), and 3968 
(solid) Earth days. All simulations adopt solar abudances for the opacities. 
Color scheme is as in Fig. [8] 

respectively, and the phase offset is ^ 57°. As can be seen 
in Fig. |7J much of this escaping radiation lies in the 3-5 /xm 
range. 

Our basic results are insensitive to the model resolution and 
integration time, as illustrated in Fig. [12] The top panel shows 
light curves at Spitzer bandpasses for solar-metallicity HD 
189733b simulations performed at horizontal resolutions of 
C64 (solid), C32 (dashed), and C16 (dotted) (approximately 
equivalent to global resolutions of 256 x 128, 128 x 64, and 
64 X 32 in longitude and latitude, respectively). All are at an 
integration time of 301 Earth days and have 40 vertical layers. 
In the bottom panel are light curves from the low-horizontal- 
resolution (C16) case at integration times of 104, 995, and 
3968 Earth days. The simulated light curves are very similar 
in all these cases. This indicates that our model resolutions 
and integration times are sufficient to capture the dynamics. 

Our simulations are consistent with recent upper limits 
on the temporal variability of HD 189733b reported by 
lAgol et all (120081) . They analyzed five Spitzer IRAC 8-^m 
secondary-eclipse observations scattered over ~ 370 Earth 
days and found that the variability in secondary-eclipse depth 
is less than 10% around their best-fit mean value of 0.00347. 
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Fig. 13. — Secondary-eclipse depths versus time calculated from our 
synchronously rotating HD 189733b simulation with solar opacity. Curves 
show results at the 3.6, 4.5, and S-iim Spitzer IRAC bands in black, red, and 
blue, respectively, each normalized to its mean. Our predicted v ariability is 
~ 1-2%, consistent with the observational 8-/xm upper limit from lAgol et alj 
<200a) . 

To compare with this observation. Fig. [13] shows secondary- 
eclipse depths in the IRAC 3.6, 4.5, and 8-//m bands ver- 
sus time, calculated from our HD 189733b simulation (as in 
Figs.[8]and[T2]l by properly integrating the simulated infrared 
spectra over the Spitzer bandpasses. Over a period of sev- 
eral hundred days, our simulated variability is ^ 1% at 4.5 
and 8 im\ and ^ 1.5% at 3.6 //m, fully consistent with the 
lAgol etalJ ('2Q08) upper limit. Interestingly, much of the pre- 
dicted variability involves a coherent oscillation with a period 
of approximately 43 days, although a weak downward trend 
(corresponding to a mean decrease in eclipse depth of about 
0.5%) is also present. The predominance of a single oscilla- 
tion period suggests the presence of a global sloshing mode 
with a period of weeks; future work will be required to inves- 
tigate this phenomenon and the extent to which it may depend 
on parameters such as the atmospheric vertical structure. 

4. NON-SYNCHRONOUSLY ROTATING HD 189733B 

Several authors have suggested that hot Jupiters in near- 
circular orbits will synchronously rotate because their ex- 
pected spindown times are ~ 1 0^ ye ars for a Jupiter-like tidal 
Q of 10^ (Guillotetal. 19961: iLubow et al. 1997). While 
this is plausible, no observations currently exist to con- 
strain the rotation period of any hot Jupiter. Moreover, 
IShowman & GuilioB(l2002l) argued that dynamical torques be- 
tween the atmosphere and interior could lead to equilibrium 
rotation rates that deviate modestly from synchronous. There- 
fore, we ran several simulations of HD 189733b to investigate 
the role that non-synchronous rotation has on the dynamics 
and resulting observables. 

Figure[T4]shows the zonal-mean zonal winds versus latitude 
and pressure for cases with rotation rates that are half, 1.5, 
and twice the synchronous values {top, middle, and bottom 
panels, respectively). These cases all assume solar metallicity 
and can be directly compared to their synchronously rotating 
counterpart in Fig. |5] As can be seen, the rotation rate has a 
significant effect on the mean jet structure. The slowly rotat- 
ing case {top) qualitatively resembles the synchronous case, 
with an eastward equatorial jet and westward high-latitude 
flows, except that the equatorial jet is narrower and faster. At 
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Fig. 14. — Effect of non-synchi'onous rotation on jet structure. Shows 
zonal-mean zonal wind versus latitude and pressure for HD 189733b cases 
with half (top), 1.5 times (middle), and twice (bottom) the synchronous ro- 
tation rate. Opacities use solar abundances. Scalebar gives speeds in ms"'. 
All simulations have a horizontal resolution of C32 with 40 layers and can be 
directly compared with Fig.|5] 



rotation rates faster than synchronous, however, the equato- 
rial jet weakens and strong eastward polar jets develop — 
a phenomena unseen in our synchronous cases. At rotation 
rates that are double the synchonous value {bottom), the mean 
speed of the polar jets exceeds that of the equatorial jet. Inter- 
estingly, these faster rotation rates are also accompanied by 
a general slowdown in the wind speeds; at double the syn- 
chronous rotation rate, the peak zonal-mean zonal wind speed 
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is 1.6kms~', half the value in the sychronous case. 

Two effects can influence the jet structure in these nonsyn- 
chronous cases. First, in the nonsynchronous cases, the il- 
lumination sweeps in longitude as seen in the planet's rotat- 
ing reference frame — in contrast to the synchronous case, 
where the illuminated region is locked to a specific range of 
longitudes. In all cases (synchronous or not) the jets are pre- 
sumably accelerated by large-scale eddies generated by the 
day-night heating gradient. This must particularly be true of 
the superrotating equatorial jet, which corresponds to a lo- 
cal maximum of angular momentum per mass and can only 
be formed by up-gradient momentum transport by eddies. 
Sweeping the heating pattern in longitude can potentially 
change the character of these jet-driving eddies and hence al- 
ter the resulting jet structure. This may explain, for example, 
the trend of decreasing jet speed with increasing rotation rate 
(Fig.O. Interestingly, the peak equatorial jet speed in an in- 
ertial reference frame — namely, the sum of the peak jet speed 
in Fig. [14] and the planetary rotation rate of fla at the equator 
— is nearly constant (to within ^ 15%) in all the cases shown 
inFigs.l5]and[T4l 

Second, changing the rotation rate changes the strength of 
the Coriolis force and, particularly, the strength of the so- 
called "/3 effect," where jS = 2^1 cos cj)/ a is the northward gra- 
dient of the Coriolis parameter For flows driven primarily by 
creation of turbulence at small scales, turbulence theory pre- 
dicts that jet widths in la titude scale as ^ 7r(t///3)^/^, where 
U is a mean wind speed ( Rhinesll975l : IVasavada & ShowmanI 
l2005h . This "Rhines" scaling suggests that slower (faster) ro- 
tation rates would lead to wider (narrower) jets. In simula- 
tions of synchronously rotating hot Jupiters. IShow man et alj 
(|2008a) showed that faster rotation indeed leads to narrower 
jets and vice versa — but the dependence of jet width on rota- 
tion rate is weaker than the inverse-square-root Rhines pre- 
diction. Presumably, this deviation occurs because, in the 
hot Jupiter case, the jets are driven not at small scales but 
by eddies induced by the day-night heating gradient, which 
have length scales comparable to the jet widths. The differ- 
ences between our nonsychronous ca ses (Fig.fT?]) and flie syn- 
chronous Q parameter variations in IShowman et alj (l2008ah 
suggest that the trends in Fig. [14] cannot be explained solely 
by the jS effect but also depend on changes in the eddy behav- 
ior as discussed above. We defer a detailed diagnosis of the 
dynamics of this process to a future study. 

Figure[T5]demonstrates that nonsynchronous rotation has a 
significant effect on the light curves. Interestingly, the phase 
offsets are smaller for the slowly rotating case (il half the 
synchronous value; top panel) and larger for the rapidly ro- 
tating cases (O 1.5 and two times the synchronous value; 
bottom panel). This behavior makes sense physically. As 
viewed in the synchronously rotating reference frame, the 
nonsynchronous rotation is equivalent to a westward motion 
of the planetary interior in slowly rotating case but an east- 
ward motion of the planetary interior in the rapidly rotat- 
ing cases. This rotation combines with the winds shown in 
Fig. I14f ' to produce (as seen in the synchronous frame) a 
strong, latitude-averaged eastward flow in the rapidly rotat- 
ing cases; in contrast, in the slowly rotating cases, the midlat- 

' ' The winds shown in Fig. 1141 are those in the planet's nonsynchronous 
reference frame; to represent the winds in the synchronously rotating refer- 
ence frame, one must add (fi — f!syn)acos(/), where f2 is the actual rotation 
rate, fisyn is the rotation rate of the synchronous reference frame, a is plane- 
tary radius, and </< is latitude. 



itudes flow strongly westward, although the equatorial jet is 
still eastward. The upshot is that the hottest regions are shifted 
eastward from the substellar point by a greater distance in the 
rapidly rotating cases than in the slowly rotating case. This 
leads to a small phase shift in the slowly rotating case and a 
large phase shift in the rapidly rotating cases (Fig.lTsll. 

Interestingly, the shape of the lightcurve is also affected — 
the flux peak is narrower in the slowly rotating case than in the 
rapidly rotating case. This seems to occur because the dayside 
hot regions span a smaller range of longitudes in the slowly 
rotating case than in the rapidly rotating cases. 

The theoretical 8-/im light curve in the slowly rotating case 
actually provides quite a good match to the observed 8-/im 
light curve, reproducing both the shape and phase offset of 
the flux peak. Nevertheless, we still do not reproduce the flux 
minimum seen after transit. The rapidly rotating cases provide 
a significantly poorer match at 8 /im, with a flux peak that is 
too broad and a phase shift that is too large. The situation 
is more ambiguous at 24 /im; both cases provide reasonable 
though not perfect matches to the observed light curve. We 
re-emphasize that these are not fits to data but rather are the- 
oretical predictions. Moreover, although it could be tempting 
to infer a planetary rotation rate from these comparisons, this 
is premature, as other factors not explored here (e.g., disequi- 
librium chemistry) could also have large effects on the light 
curves. 

5. SYNCHRONOUSLY ROTATING HD 209458B 
5.1. Circulation regime, spectra, and light curves 

Spitzer secondary-eclipse photometry shows that the 4.5 
and 5.8-/im bri ghtness temperature s of HD 209458b reach 
1700-1900K (.Knutson et al.l l2008h . These values signif- 
icantly exceed both the planet's effective temperature and 
the inferred brightness temperatures at the 3.6, 8, and 24- 
/im bands, wh ich are clos er to ^ 1500K . Given these facts, 
iKnutson et all (|2008) and iBurrows et all (l2007l |2008) sug- 
gested that the atmosphere of HD 209458b contains a thermal 
inversion (hot stratosphere) on the dayside, with temperatures 
reaching ^ 2000 K or more. 

Lon g be fore these Spitzer data were available. lHubenv et al.l 
( |2003|) and lFortnev et al.l (l2006b ) showed that such inversions 
result naturally from gaseous TiO and VO, which are ex- 
tremely opaque in the visible wavelength range and cause ab- 
sorption of starlight at substantially lower pressures (■^mbar) 
than would occur in the absence of these species. Unde- 
termin ed photochemical products are another p ossible ab- 
sorber ( IBurrows et al.ll20()7l 120081) . IZahnle et al.l (12009), for 
example, have suggested that disequiUbrium sulfur species 
produced from photochemical destruction of H2S could pro- 
vide sufficient short-wavelength absorption to account for the 
stratospheres. 

However, to date, radiative-equilibrium models of HD 
209458b ha ve been unable to re produce the observed Spitzer 
photometry dFortnev et al.ll2008 ): while they successfully re- 
produce the high 4.5 and 5.8-/j,m fluxes, they overpredict the 
3.6, 8, and 24 /im fluxes. This suggests that the stratosphere 
predicted in these radiative-equili brium models is too hot 
and/or too broad in vertical extent. IBurrows et al.l (120071) ob- 
tained better agreement, especially at 3.6 /xm, by including an 
ad hoc heat sink to mimic the effect of day-to-night heat trans- 
port by the atmospheric circulation. By confining this sink to 
pressures between 0.01 and 0.1 bars. Burrows et al.l (l2007h 
were able to reproduce the low 3.6-/xm flux while keeping 
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Fig. 15. — Light curves for ou r no n-sychronous, solar-metallicity HD 
189733b simulations shown in Fig. 1 141 Top panel shows case with half the 
synchronous rotation rate (squares, asterisks, and triangles denote different 
times). Bottom panel shows cases with 1.5 times the sychronous rotation rate 
(triangles and asterisks denoting different times) and twice the synchronous 
rotation rate (diamonds). 



the 4.5 and 5.8-/im fluxes high. However, as IShowman et al.l 
(l2008a) pointed out, day-night heat transpo rt is unlikely to 
be confi ned to a narrow pressure range; in IShowman et al.l 
(l2008ai) 's 3D circulation models, for example, the dayside 
"heat sink" (expressed in Ks"' ) caused by dynamics increases 
monotonically with decreasing pressure from ^10 bars to the 
top of their model at 0.001 bar (their Fig. 10). Thus, an im- 
portant question is whether a 3D dynamical model, coupled 
to realistic radiative transfer and including TiOA'O opacity, 
can match the observed secondary-eclipse spectrum of HD 
209458b. 

Light curves are also of interes t. From sever a l brief Spitzer 
observations at different phases, ICowan et al.l (l2007h placed 
a 2-(T upper limit of 0.0015 on the peak-to-peak phase varia- 
tion in the planet/star flux ratio at 8 /im. This suggests that 
the difference in the planet's day and night 8-/im bright- 
ness temperatures is less than a few hundred K. In contrast, 
iFortney et al.l (l2008h suggested that, as a pM-class planet, HD 
209458b should exhibit much greater phase variations in the 
mid-infrared than planets lacking atmospheric TiO/VO. Phys- 
ically, the idea is that the presence of TiO and VO moves 
the photospheres upward to low pressure, where the radia- 
tive time constants are short and the air can experience rapid 



dayside heating and nightside cooling. In contrast, plan- 
ets without atmospheric TiO/VO would have deeper photo- 
spheres, where the radiative time constants are longer than 
typical advection ti mes, thus yiel ding more modest day-night 
phase variations. IFortney et al.l (|2008) also suggested that 
phase offsets of the hottest/coldest regions from the substel- 
lar/antistellar points would be smaller for pM-class planets 
than for pL-class planets. We here wish to test these ideas 
by determining the circulation patterns and lightcurves for 
a model of HD 209458b with TiO/VO opacity. Note that 
even if hot-Jupiter stratospheres turn out to result from short- 
wave absorption b y compounds other than TiO and VO (e.g., 
IZahnle et alfmO^ . our simulations will give a qualitative pic- 
ture of how dynamics responds to the presence of a compound 
that absorbs strongly in the visible. 

Thus, we performed several simulations of HD 209458b in- 
cluding TiO and VO opacity in chemical equilibrium. When 
temperatures are too cold for TiO and VO to exist in the 
gas phase, they are absent, but when temperatures are warm 
enough, they are included. 

Our HD 209458b simulations develop vigorous circulations 
that, as expected, include a hot dayside stratosphere. This is 
illustrated in Fig. [16] which shows the temperature and hor- 
izontal wind patterns at pressures of 1 bar, 30 mbar, 1 mbar, 
and 0.1 mbar {bottom to top) for our case with solar abun- 
dances. At deep levels (1 bar and 30 mbar), the circulation 
qualitatively resembles that of HD 1 89733b, with an eastward 
equatorial jet, westward high-latitude flows, and a hot region 
shifted east of the substellar point (compare bottom two pan- 
els of Figs.l4land[T6]l. Nevertheless, even at these deep levels, 
localized regions attain temperatures sufficient for TiO/VO to 
exist in the gas phase (- 1700-1900 K). 

By pressures of 1 mbar, the picture changes drastically: 
a "baby" stratosphere with a radius of 50° in longitude and 
latitude has developed, with temperatures reaching 2000 K 
(Fig. [16] second panel). This stratosphere is approximately 
centered on the substellar point, with an eastward offset of 
only 10° longitude. Its spatial confinement results from 
the fact that only air within ~ 50° of the substellar point re- 
ceives sufficient irradiation to achieve temperatures necessary 
for gas-phase TiO and VO to exist. Air > 60° from the sub- 
stellar point, although still on the dayside, has temperatures 
too low for gas-phase TiO/VO, and without the benefit of the 
extra opacity supplied by these species, this air remains sub- 
stantially cooler (< 1400 K). This nonlinear feedback between 
temperature and opacity leads to a remarkably sharp tempera- 
ture gradient, with temperatures dropping from 1800 to 1400 
K as one moves from 50° to 60° angular distance from the 
substellar point. 

At still lower pressures, an even greater fraction of the 
starlight is available to cause heating, and the stratosphere be- 
comes horizontally larger until it covers most of the dayside 
at pressures < 0.1 mbar (Fig. [16] top panel). Although the 
equatorial winds at these low pressures involve a significant 
eastward flow at most longitudes, the mid- and high-latitude 
winds involve a simpler motion that, to zeroth order, moves air 
away from high-temperature regions toward low-temperature 
regions. At low pressure, the temperature pattern is relatively 
symmetric about the substellar and a ntistellar points, as sug- 
gested observationall y for Ups And b ([Harrington et al.l2006l) 
and HD 179949 (Co wan et al.ll2007h . 

A contributing factor to the widening of the stratosphere 
with altitude is that stellar radiation is absorbed at lower pres- 
sure near the limb than at the substellar point, which results di- 
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Fig. 16. — Temperature (colorscale, in K) and winds (aiTows) for nomi- 
nal HD 209458b simulation with solar abundances including TiOA'O. Panels 
show flow at 0.1 mbar (top), 1 mbar {second panel), 30 mbar (third panel), 
and 1 bar (bottom panel). Horizontal resolution is C32 (roughly equivalent to 
a resolution of 128 X 64 on a longitude/latitude grid) with 53 vertical layers. 
Substellar point is at longitude, latitude (0° , 0° ). Note the development of the 
dayside stratosphere. 
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Fig. 17. — Chemical-equilibiium abundances of several molecules versus 
pressure at the substellar point (solid) and antistellar point (dashed) from our 
three-dimensional HD 209458b case adopting solar abundances. 



rectly from the longer atmospheric path length to stellar irra- 
diation at the limb (ji 0) than at the substellar point (/i = 1). 
Thus, near the limb the heating peaks at lower pressure, con- 
fining the stratosphere to low pressures. Near the substellar 
point, the heating peaks deeper, allowing the stratosphere to 
extend more deeply in this localized region. 

Our simulated stratospheres result directly from the exis- 
tence of gaseous TiO and VO on the dayside in these simula- 
tions. Simulations of HD 209458b performed without gaseous 
TiO/VO lack stratospheres and develop temperature and wind 
patterns resembling a hotter version of our HD 189733b sim- 
ulations. Likewise, our previous HD 209458b simulations 
per formed with Newtonian heating/cooling (Showman et 
l200 8a) developed a weak dayside temperature inversion but 
nothing resembling the extremely hot stratosphere shown in 

To illustrate the height-dependence of the composition 
(which affects the opacity). Figure [17] depicts the chemical- 
equilibrium abundances versus pressure at the substellar point 
(solid) and antistellar point (dashed) for our nominal HD 
209458b case. At low pressure, TiO is abundant on the day- 
side yet depleted on the nightside. 

Importantly, our simulated stratospheres do not extend fully 
to the terminators (which are at longitudes ±90° in Fig.fTSI). 
Generally, the terminators themselves have temperatures of 
~ 1300K or less. The net heating there is simply too low to 
allow temperatures above the TiOA'^O condensation curves. 
Thus, while much of the dayside could have abundant TiO and 
VO, our simulations suggest that these species will be largely 
absent (or at least depleted) on the limbs as seen during pri- 
mary transit. This will have important impli cations for inter- 
preting transit spectra of HD 209458b (e.g. Sing et al.ll2008h 
and other pM-class planets. Interestingly, note that our simu- 
lated stratosphere approaches the teminator to the east of the 
substellar point more closely than it approaches the teminator 
to the west of the substellar point — the result of thermal ad- 
vection due to the eastward equatorial jet. This result suggests 
that, during transit, the planet's leading limb will be cooler 
and more depleted in TiOA' O than the trailing limb. 

The diversity of temperature-pressure profiles in our solar 
model is shown in Fig. [18] Near the bottom is the quasi- 
isothermal region from ^ 3-100 bars. Above that, from pres- 
sures of 1 bar to 30 mbar, the temperature decreases with 
altitude on both the dayside and the nightside. At pressures 
less than ~ lOmbar, the dayside stratosphere is evident, with 
temperatures exceeding 2000 K. 
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Fig. 18. — A selection of temperature-pressure profiles for our nomi- 
nal, solar-abundance HD 20 945 8b simulation including TiO and VO opacity 
(same simulation as in Fig. |16) . Red (blue) profiles ai'e equatorward (pole- 
ward) of 30° latitude. Note the formation of the stratosphere at pressures less 
than ~ 30mbar. 



Figure[T9]shows spectra at six orbital phases. Within 60° 
orbital phase of transit {magenta, black, and red curves), 
molecular bands are seen in absorption because the tempera- 
ture decreases with increasing altitude on the nightside. Near 
secondary eclipse, however, when the dayside faces Earth, 
the features flip into emission (green, dark blue, and light 
blue curves), which results directly from the dayside tem- 
perature inversion associated with the stratosphere. This dif- 
fers from our HD 189733b simulations, where molecular fea- 
tures are seen in absorption at all orbital phases, including 
the dayside (Fig. |7ll. Nevertheless, as the effective temper- 
ature of this HD 209458b simulation makes clear (Fig. [TTT i. 
the dayside radiation arrives primarily from the bottom of the 
stratosphere where mean temperatures are a modest ^ 1500K 
(rather than from higher-altitude regions where temperatures 
exceed 2000 K). 

Figure |20] compares our simulated HD 209458b dayside 
planet/star flu x ratio spectrum to \hs Spitze r sec ondary-echpse 
photometry of iKnutson et al.l (l2008l l2009l) and lDeming et al.1 
OP05), including a tentative revision to the 24 /xm eclipse 
depth by Deming (personal communication). We match the 
secondary-eclipse depths at 3.6 and 8/im. However, de- 
spite the existence of a stratosphere in our simulations, our 
solar-opacity case underpredicts the eclipse depths at 4.5 and 
5.8/im by nearly 50% (roughly 3 and 23-a, respectively, at 
these two b ands). We also miss the 24-fim eclipse depth of 
iDeming eT al. (2005) (the lower point in Fig. |20] | by roughly 
2.5-(j, although we fall within 2-cr if the tentative revision to 
this point (upper point in Fig. |20] | turns out to be more appro- 
priate. 

These discrepancies suggest that our simulated stratosphere 
does not have the correct properties (e.g., temperature or 
altitude range). A comparison of our temperat ure profiles 
(Fig.fTsTl to the radiative-equilibrium profiles of Fortney et alJ 
(12008, their Fig. 12), who match the 4.5 and 5.8-/zm points 
but overpredict the others, shows that the hottest of our strato- 
spheric profiles qualitatively resembles theirs. However, be- 
cause our dayside includes a range of profiles ranging from 
hot to cold, our dayside average profile is colder than theirs. 
This may help explain the fact that our predicted dayside 
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Fig . 19. — Emergent flux density (ergs s cm Hz" ) from our nominal, 
solar-abundance simulation of HD 209458b, including TiO and VO opacity, 
at six orbital phases. Black, nightside, as seen during transit; red, 60° after 
transit; green, 120° after transit; dark blue, dayside, as seen during secondary 
eclipse; light blue, 60° after secoSdary eclipse; and magenta, 120° after sJe- 
ondaiy eclipse. The key in the top right comer is color-coded with the spectra 
to illustrate the sequence. Thin dotted black lines at the bottom of the figure 
show normalized Spitzer bandpasses and the letters at the top show loc ations 
of the H, K, L, and M bands. This is the same simulation as in Fig. 1161 



fluxes are systematically lower than those of iFortnev et alj 
( |2008|) . Presumably, this difference results from the vigor- 
ous dayside-to-nightside transport of thermal energy by the 
atmospheric circulation in our case. 

When a spectrum exhibits greatly differing brightness 
temperatures at different wavelengths (as is true for HD 
209458b), a common explanation is that the different wave- 
lengths sense different pressure levels (because of the 
wavelength-dependent opacities). In the presence of a ver- 
tically varying temperature this can produce a spectrum with 
wavelength-varying brightness temperatures. 

In this context, a fundamental stumbling block to simulta- 
neously explaining the five Spitzer secondary-eclipse depths 
is that the range of pressures that contribute photons to the 3.6, 
4.5, 5.6, 8, and 24-/im bands are all very similar — at least for 
the radiative-transfer model, chemical composition, and opac- 
ities adopted here. Contribution functions calculated for each 
of these Spitzer bands (see Fig.|9l right panel) peak between 
3 and 10 mbar, and they all have very broad tails extending 
to ^ 0.1 mbar on the low-pressure side and ^ lOOmbars on 
the high-pressure side.'^ This overlap in pressure means the 
brightness temperatures in all these bands tend to be similar. 
Thus, it is difficult to produce a high brightness temperature in 
some bands (such as at 4.5 and 5.8 pm) while maintaining low 
brightness temperature in other bands (such as at 3.6, 8, and 
24 pm) — as apparently required by the data. This problem is 
not confined to the pr esent study; it helps explain the difficulty 
IFortnev et alj (12008 ) had i n explaining all the observations. 
Even lBurrows et aL (l2007l) . who had the flexibility of several 
free parameters governing an assumed stratospheric absorber 
and the magnitude and pressure range of the parameterized 
dayside heat sink, were unable to match the 5.8/8-/im flux 
ratio. Separating the pressure ranges of the contribution func- 

These pressure ranges ai'e specific to the solar model including TiO 
and VO. Contribution functions calculated without TiOA^O for HD 189733b 
using the same radiative-transfer code peak at substantially deeper pressures 
of 30-lOOmbai-. SeeFig.|9] 
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Fig. 20. — Planet-to-star flux ratio vs. wavelength for our HD 209458b 
simulation at the time immediately before/after secondary eclipse {solid 
cun'e). The simulation adopts solar abundances including TiOA'O opac- 
ity. Points show measure d secondary-eclipse depths tPeming et aTl |200^ 
IKnutson et aU2008ll200^ . 

tions would require differential changes to the opacities at the 
wavelengths of Sp/fzer bandpasses. Potentially, alternate (i.e., 
disequilibrium) chemical compositions could sufficiently af- 
fect the opacities to resolve this conundrum; investigating this 
possibility will require future work. Disequilibrium chemistry 
has been invoked to explain infrare d spectra of L d warfs with 
similar effective temperatures (Leg gett et al.ll2007l) . 

We now turn to infrared light curves. Figure |2T| shows 
light curves in Spitzer bandpasses calculated for our HD 
209458b case with solar abundances. As was the situation 
with HD 189733b, our current HD 209458b light curves 
exhibit muted phase variations relative to those obtained 
from our earlier simulat i ons using Newtonian heating/cooling 
dShowman et al.ll2008at iFortnev et aT]|2006al) . In the current 
simulations, the ratio of maximum-to-minimum planet/star 
flux ratio (at a given wavelength) is nearly 2 at all four IRAC 
bands and 1 .6-1 .7 at the 16 /im IRS and 24 /im MIPS bands. 
The difference between maximum and minimum values in- 
creases with wavelength and is 0.0014 at 8 fim. Phase offsets 
range from ^ 20° at the longer wavelengths to almost 40° at 
3.6 fim. These values are significantly nonzero but are smaller 
than their counterparts for our solar-abundance HD 189733b 
simulations (Fig. |8l top panel), consis tent with the general 
trend suggested bv lFortnev etal](l2008l) . 

Our muted phase variations are consistent with the obser- 
vational constraint of Cowan et al. (2007), who found a 2-a 
upper Hmit of 0.0015 on the phase variation of HD 209458b 
at 8 fim. Fig. |2T]also re-iterates our agreement with the 3.6 
and 8 /im secondary-eclipse photometry and our discrepancy 
at 4.5 and 5.8 /im. Continuous light curves over half an orbit 
have recently been obtained at 8 and 24 /im and are possible 
at 3.6 and 4.5 /im with warm Spitzer, which would provide a 
test of our models and insights on how to explain the discrep- 
ancies noted above. 

The muted phase variations in Fig. |2T|are particularly in- 
triguing given that we do have a dayside stratosphere and (as 
a result) enormous day-night temperature variations reaching 
~ 1500K at low pressure. There are two possible reasons 
to reconcile why these large day-night temperature variations 



do not translate into large day-night flux variations. First, 
much of the radiation that contributes to the dayside fluxes 
emanates from altitudes near the bottom of the stratosphere, 
where day-night temperature differences are modest. This is 
exacerbated by the fact that the stratosphere only covers part 
of the dayside. Thus, a substantial fraction of the dayside flux 
does not emanate from the stratosphere but from deeper layers 
and surrounding regions that are cooler. This lowers the day- 
side flux compared to an idealized case where all the dayside 
mid-infrared radiation comes from the stratosphere. Second, 
the cooler temperatures on the nightside means that the mid- 
IR photospheres on the nightside are deeper in pressure — by 
roughly an order-of-magnitude — than those within the day- 
side stratosphere. Because temperatures increase with pres- 
sure on the nightside, this elevates the nightside flux relative to 
an idealized case where the photospheres remain at the same 
(low) pressure everywhere. Both these effects act to mute the 
phase variations in the mid-lR, at least in the Spitzer band- 
passes. 

Note, however, that the expected IR phase variations can 
depend sensitively on wavelength. As can be seen in Fig. [19] 
for example, our HD 209458b model predicts that the peak 
phase variations reach factors of 4 in specific wavelength 
bands near 1.5, 1.8, and 2.8 /im, with relatively smaller phase 
variations (factors of ^ 2 or less) at intermediate wavelengths. 
Since both large and small day-night phase variations can oc- 
cur on a single object (depending on wavelength), caution is 
needed when attempting to link the amplitude of phase varia- 
tion to the efficiency of day-night heat transport — especially 
if one only possesses light curves at only one or two wave- 
lengths. Light curves obtained in numerous isolated bands 
across this region (or low-to-moderate resolution spectra at 
different orbital phases) by a future space mission could shed 
considerable insight into the atmospheric composition and 3D 
temperature structure. 

5.2. Discussion 

Our results p oint toward possible refinements of the sce- 
nario outlined in lFortnev et al.l ( l2008l) . They anticipated large 
day-night flux differences for pM-class planets, including HD 
209458b. This contrasts with the modest phase variations we 
find here. However, the fact that we underpredict the day- 
side 4.5 and 5.8-/im fluxes for HD 209458b suggests that in 
reality (as opposed to in our simulations) the day-night flux 
variations will be large at these two wavelengths. If similar 
flux ratios are found for other planets with da y-side tempera- 
ture i nversions (XO-lb m ay point to th i s: see [ Machalek et^ 
l2008h . a revision to the iFortney et al.l (l2008h scenario may 
be in order, wherein HD 209458b and other transitional pM- 
class planets will have large phase variations at 4.5 and 5.8 /im 
but modest phase variations at the other Spitzer bandpasses. 
K band, sensing even deeper than 100 mbar, should also 
exhibit only modest phase variations. This scenario recon- 
ciles the secondary-eclipse spectrum of HD 209458b with 
the observational u pper limit on the phase variation at 8 fim 
dCowanet al.ll2007h . 

But what of hotter planets? F ortnev et al.l (|2008) posi- 
tioned HD 209458b near the pL/pM boundary and acknowl- 
edged that this transition could be indistinct because of the 
extended temperature range over which TiO and VO con- 
dense (Lodders 2002). Planets more strongly illuminated 
than HD 209458b should have hotter stratospheres that cover 
more of the dayside, and it remains possible that such "very 
hot Jupiters" could indeed have large day-night phase varia- 
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Fig. 21. — Light curves versus orbital phase calculated in Spitzer band- 
passes for our HD 209458b case with solai' abundances. Within each panel, 
moving from bottom to top, the light curves are for wavelengths 3.6 /^m 
(black), 4.5 /im (red), 5.8 ^m (green), 8/im (dark blue), 16/im (light blue), 
and 24 /im (magenta), resp ectively. Overplotted are the Spitzer secondary- 
eclipse measurements from lKnutson et all j2008l) in large diamonds, color- 
coded to the light curves. 

tions i n most Spitzer channels as predicted by iFortnev et alj 
( 20081). The strong o bservational evide nce that ups And b 
faarrington et al.l200 6) and HD 179949b (ICowan et al.l2007h 
show large day-night flux contrast at 24 and 8/im, respec- 
tively, may indicate such a change in temperature structure at 
the higher incident fluxes 35% higher than HD209458b) 
that the se planets receive. Although neither of these planets 
echpse. lB^nes et al.l(l2008 h have suggested that HD 179949b 
has a dayside temperature inversion. It will be interesting 
to perform 3D simulations of these and hotter planets, par- 
ticularly ones that go through secondary eclipse. Obtaining 
light curves for very hot eclipsing systems, such as TrES-4 
and HAT-P-7b, should be important goals for warm Spitzer or 
JWST. 

For all of these systems, whether large day-night flux con- 
trasts in the IR translate into large day-night effective tem- 
perature contrasts can only be unambiguously answered by 
light curves across many wavelengths, pref erably from 2- 
5 /Ltm, around the peak in planetary flux (e.g. lBarmanll2008h . 
This makes 3.6 and 4.5 /im light curves from the Spitzer warm 
mission, and continued searches for planet flux in K band, par- 
ticularly important. 

Finally, we conclude this section with some discussion of 
the likelihood that gaseous Ti O and VO ca n actua lly exist on 
the dayside. As discussed in IFortnev et al.l (l2008l) . we calcu- 
late opacity assuming local chemical equilibrium at the given 
p and T . Assuming TiO and VO are included in the database 
in the first place, this means that gaseous TiO/VO opacity are 
included if the local temperatures are hot enough — ignoring 
the possibility of any "cold trap" effect that could globally de- 
plete the abundance of these species. One possible cold trap, 
commonly discussed in the context of ID models, is at pres- 
sures of tens to hundreds of bars where the temperature at 
the bottom of the near-isothermal radiative zone is expected 
to be cooler than the TiO/VO condensation temperatures for 
multi- Gyr-old hot Jupiters farther than <y 0. 04 AU from their 
stars ( Hubeny et al.ll2003HFortney et"aD '2008). Nevertheless, 
the existence of stratospheres on some planets argues that this 
does not occur on HD 209458b and other pM-class planets. 



However, another possible cold trap — not considered by 
previous ID models — is the presence of the large day-night 
temperature difference in the observable atmosphere. As a 
parcel of hot dayside air flows onto the nightside, its temper- 
ature plummets and gaseous TiO and VO condense. If these 
Ti- and V-bearing particles settle out before the air parcel re- 
turns to the dayside, then the atmosphere could become de- 
pleted in TiO and VO even if no cold trap exists at deeper 
levels of tens to hundreds of bars. On the other hand, if the 
particles are small, their settling speeds will be slow and they 
will remain in the air parcel when it returns to the dayside 
days later. In this case, the increased temperatures will al- 
low these particles to sublimate, resupplying gaseous TiO and 
VO to the atmosphere. Thus, our simulations with TiO/VO 
are only self-consistent if the TiO/VO particles cannot settle 
out on the nightside. This effectively means that the particles 
must remain small. 

We can quantify the maximum particle sizes as follows. At 
1 mbar, characteristic vertical velocities in our HD 209458b 
simulation are ^ 30ms~'. For particles to remain suspended, 
the settling velocities must be less than these values. Given 
the Stokes flow speed modifie d for gas-kinetic effects (see 
e.g. lAckerman & Marleyl 12001 L Eq. Bl), and air viscosities 
relevant to hydrogen at ~ 1500K, this implies particles radii 
less than ^ 30 //m. At lower pressure, the gas-kinetic effects 
become stronger, leading to greater settling velocities; at 0.1 
mbar, for example, the particles must be less than ^ 7-10/im 
in radius to remain suspended. (Calculation of Reynolds num- 
bers for the falling particles shows that they are less than one, 
implying that the Stokes relation is valid and that turbulent 
modifications to the fall velocity need not be considered.) It 
seems plausible that the actual particles sizes are smaller than 
these values, but detailed microphysical calculations will be 
needed to explore this further. A complicating factor is that 
not only TiO/VO but also silicates will condense on the night- 
side at low pressure. 

These estimates are c rudely consist ent with subsequent cal- 
culations performed bv lSpiegel et al.l ( 12009.) . who parameter- 
ized the dynamics in ID using eddy diffusion and estimated 
the eddy-diffusion coefficients needed to keep TiO particles 
of various sizes lofted. If one translates our vertical velocities 
into eddy diffusivities by multiplying the vertical velocities by 
the pressure scale height, then the implied eddy diffusivities 
at 1 mbar ai-e 10"cm^sec"'. In their HD 209458b model, 
at 1 mbar, these eddy diffusivities are sufficient to loft parti- 
cles smaller than several /im. Thus, their calculations likewise 
suggest that TiO stratospheres may be viable if the TiO parti- 
cle sizes are less than several /im. 

6. CONCLUSIONS 

We presented global, three-dimensional numerical simula- 
tions of the atmospheric circulation of HD 189733b and HD 
209458b that couple the dynamics to a realistic representation 
of cloud-free non-gray radiative transfer This new model, 
which we dub SPARC/MlTgcm, is the first 3D dynamical 
model for any giant planet — including those in our Solar Sys- 
tem — to incorporate nongray radiative transfer Our model 
adopts the MlTgcm for the dynamics and an optimized ver- 
sion of the radiative model of McKay, Marley, Fortney, and 
collaborators for the radiative transfer Opacities are calcu- 
lated assuming solar composition (or some multiple thereof) 
with equilibrium chemistry that accounts fo r rainout. Like 
earlier work with simplified forcing ( Showman et al.l|2008al; 
ICooper & ShowmanI l200l IShowman & GuillotI |2002|) . our 
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simulations develop a broad eastward equatorial jet, mean 
westward flow at high latitudes, and substantial flow over the 
poles at low pressure. The jet structure depends significantly 
on longitude at pressures < lOOmbar. 

For HD 189733b, our simulations that exclude TiO and 
VO opacity explain the broad features of the ob served 8 
and 24-/im light curves dKnutson et al.ll2007l l2009h . includ- 
ing the modest day-night flux variation and the fact that the 
planet/star flux ratio peaks before the secondary eclipse. In 
our simulations, the offset results from the eastward displace- 
ment of the hot regions from the substellar point. On the other 
hand, we do not fit the flux minimum seen after transit in the 
8-/im light curve (iKnutson et al.ll2007h . Our simulations also 
provide reasonable matches to the Spitze r secondary-eclipse 
depths at 4.5, 5.8, 8, 16 , and 24 /im dCharbonneau et alj 
l2008t iDeming et al.l 12006) and the groundbased upper Umit 
at 2.2 ^m from lBarnes et al.l (l2007l) . The temporal variability 
in these simulations is modest — of order 1 % — and is fully 
consistent with th e upper limit on temporal variability from 
lAgoletalJdlOOSh . 

The primary HD 189733b observation where we 
fare poorly is the 3.6j!i m secondary-eclipse depth from 
ICharbonneau et all (l2008h . which we underpredict by about 
a factor of two. Because the 3.6 /im channel is expected to 
sense down to ^ 0.1-1 bar on this planet (Fig.|9] left panel), 
this suggests that our simulation is too cold in this region 
of the atmosphere and/or has the incorrect opacity in this 
wavelength range. Previous ID models of HD 189733b have 
suffered a similar pr oblem at this wavelength (lBarmanll2008l : 
iKnutson et alJ 12009*). as have ID models of brown dwarfs 
with similar effective temperatures (Geballe et al. 2009). A 
full-orbit light curve of HD 189733b at 3.6 /im, possible with 
warm Spitzer, would provide crucial insights to help resolve 
this problem. 

For HD 209458b, we include gaseous TiO and VO opac- 
ity to see whether it allows us to explain the inference of 
a stratosphere from 5pfeer photometry (fKnutsonet al.H20d8i : 
iBurrows et"ani2007l:lFortnev et alj|2008h . As expected, these 
simulations develop a hot (> 2000 K) day side stratosphere 
whose horizontal dimensions are small at depth but widen 
with altitude until the stratosphere covers most of the day- 
side at pressures < O.lmbar. Interestingly , both b ranches of 
the bifurcation discussed by iHubeny et al.l (l2003h for differ- 
ent planets occur here on a single pla net's dayside. Using 
the terminology of lFortnev et alj (l2008h . the substellar region 
is a pM-class planet but the terminators and nightside are a 
pL-class planet. It is thus perhaps more proper to talk about 
pM-class daysides rather than pM-class planets. 

But despite the stratosphere in our simulations of HD 
209458b, we do not reproduce current Spitzer photometry of 
this planet, which includes particularly high 1700-1900 K) 
brightness temperatures in the 4.5 and 5.8-/im channels. This 
could mean that our stratosphere has the incorrect proper- 
ties (e.g., temperature range, altitude range, and vertical ther- 



mal gradient). However, a fundamental difficulty in explain- 
ing the diverse brightness dayside temperatures (ranging from 
~ 1500-1900K in Spitzer bandpasses) is that the range of 
pressures from which emergent infrared photons originate 
are very similar for all Spitzer bandpasses, at least as calcu- 
lated by our radiative model with equilibrium chemistry. This 
means that regardless of the planet's temperature profile, the 
brightness temperatures in all these bands should be similar 
Breaking this degeneracy would require changing the opaci- 
ties so that the opacities in different Spitzer bandpasses differ 
significantly. Dropping the equilibrium chemistry assumption 
would make this possible; this provides a clue that disequi- 
librium chemistry may be important. Disequilibrium chem- 
istry ap pears to influence infrared spectra o f brown dwarfs 
(e.g. Leggett et al.l l2007l: iGeballe et al.ll2009l) . lending weight 
to this possibility. 

Our light curves of HD 209458b in Spitzer bandpasses 
exhibit modest day-night variation, and we successfully ex- 
plain the up per limit on the day-night flux contrast from 
ICowan et al.l (1200 7) at 8 pm. Our ability to meet this con- 
straint results from the fact that much of the dayside 8-/im 
radiation emanates from altitudes near the base of the strato- 
sphere, where temperatures are not too hot, while nightside 
8-/im radiation emanates from substantially deeper pressures 
where temperatures are warmer than they are aloft. This dual 
effect leads to modest day-night flux variations despite large 
day-night temperature variations (on isobars) at low pressures. 
Assuming the high inferred 4.5 and 5.%-pm. brightness tem- 
peratures indeed result from a stratosphere, we suggest that 
the real planet should exhibit large phase variations at 4.5 
and 5.8 /im yet more modest phase variations at other Spitzer 
bandpasses and K band. 

The task of developing exoplanet GCMs that couple dy- 
namics and radiative transfer has now been espoused by 
many authors, and the work presented here demonstrates 
that this approach indeed holds significant promise for ex- 
plaining atmospheric observations of hot Jupiters. Our 3D 
SPARC/MITgcm simulations, especiaUy of HD 189733b, 
show encouraging resemblance to observations of the real 
planet. While some discrepancies remain, and a wider range 
of parameters need be explored, these simulations support 
the idea that detailed model/data comparisons can eventually 
allow robust inferences about the circulation regime of hot 
Jupiters to be inferred. Given the huge range in properties of 
currently known transiting planets, with future observations 
sure to unveil additional surprises, this helps energize the ex- 
citing prospect that planetary meteorology can successfully 
be extended beyond the confines of our Solar System. 
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